songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ENH: handle uniform Su/Sp/SuSp source terms directly

Browse Source 
- avoid any operations for zero sources

- explicit sources are applied to the entire mesh can be added directly,
  without an intermediate DimensionedField

- update some legacy faMatrix/fvMatrix methods that used Istream
  instead of dictionary or dimensionSet for their parameters.
  Simplify handling of tmps.

- align faMatrix methods with the updated their fvMatrix counterparts
  (eg, DimensionedField instead of GeometricField for sources)
This commit is contained in:
Mark Olesen 2022-05-27 17:13:51 +02:00
parent 21234ae296
commit 96cc6024c0
17 changed files with 1380 additions and 917 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/faOptions/faOption/faOptionList.H
View File

@ -34,8 +34,8 @@ SourceFile
\*---------------------------------------------------------------------------*/
#ifndef faOptionList_H
#define faOptionList_H
#ifndef Foam_faOptionList_H
#define Foam_faOptionList_H
#include "faOption.H"
#include "PtrList.H"
@ -46,9 +46,11 @@ SourceFile
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Forward Declarations
namespace Foam
{
// Forward Declarations
namespace fa
{
class optionList;

4
src/finiteArea/faMatrices/faMatricesFwd.H
View File

@ -28,8 +28,8 @@ Description
\*---------------------------------------------------------------------------*/
#ifndef faMatricesFwd_H
#define faMatricesFwd_H
#ifndef Foam_faMatricesFwd_H
#define Foam_faMatricesFwd_H
#include "fieldTypes.H"

580
src/finiteArea/faMatrices/faMatrix/faMatrix.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2019-2021 OpenCFD Ltd.
Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -200,26 +200,18 @@ Foam::faMatrix<Type>::faMatrix
<< endl;
// Initialise coupling coefficients
forAll(psi.mesh().boundary(), patchI)
forAll(psi.mesh().boundary(), patchi)
{
internalCoeffs_.set
(
patchI,
new Field<Type>
(
psi.mesh().boundary()[patchI].size(),
Zero
)
patchi,
new Field<Type>(psi.mesh().boundary()[patchi].size(), Zero)
);
boundaryCoeffs_.set
(
patchI,
new Field<Type>
(
psi.mesh().boundary()[patchI].size(),
Zero
)
patchi,
new Field<Type>(psi.mesh().boundary()[patchi].size(), Zero)
);
}
@ -236,7 +228,6 @@ Foam::faMatrix<Type>::faMatrix
template<class Type>
Foam::faMatrix<Type>::faMatrix(const faMatrix<Type>& fam)
:
refCount(),
lduMatrix(fam),
psi_(fam.psi_),
dimensions_(fam.dimensions_),
@ -246,75 +237,56 @@ Foam::faMatrix<Type>::faMatrix(const faMatrix<Type>& fam)
faceFluxCorrectionPtr_(nullptr)
{
DebugInFunction
<< "copying faMatrix<Type> for field " << psi_.name()
<< endl;
<< "Copying faMatrix<Type> for field " << psi_.name() << endl;
if (fam.faceFluxCorrectionPtr_)
{
faceFluxCorrectionPtr_ = new
GeometricField<Type, faePatchField, edgeMesh>
(
*(fam.faceFluxCorrectionPtr_)
);
faceFluxCorrectionPtr_ =
new GeometricField<Type, faePatchField, edgeMesh>
(
*(fam.faceFluxCorrectionPtr_)
);
}
}
template<class Type>
Foam::faMatrix<Type>::faMatrix
(
const GeometricField<Type, faPatchField, areaMesh>& psi,
Istream& is
)
Foam::faMatrix<Type>::faMatrix(const tmp<faMatrix<Type>>& tmat)
:
lduMatrix(psi.mesh()),
psi_(psi),
dimensions_(is),
source_(is),
internalCoeffs_(psi.mesh().boundary().size()),
boundaryCoeffs_(psi.mesh().boundary().size()),
lduMatrix(tmat.constCast(), tmat.movable()),
psi_(tmat().psi_),
dimensions_(tmat().dimensions_),
source_(tmat.constCast().source_, tmat.movable()),
internalCoeffs_(tmat.constCast().internalCoeffs_, tmat.movable()),
boundaryCoeffs_(tmat.constCast().boundaryCoeffs_, tmat.movable()),
faceFluxCorrectionPtr_(nullptr)
{
DebugInFunction
<< "constructing faMatrix<Type> for field " << psi_.name()
<< endl;
<< "Copy/Move faMatrix<Type> for field " << psi_.name() << endl;
// Initialise coupling coefficients
forAll(psi.mesh().boundary(), patchI)
if (tmat().faceFluxCorrectionPtr_)
{
internalCoeffs_.set
(
patchI,
new Field<Type>
(
psi.mesh().boundary()[patchI].size(),
Zero
)
);
boundaryCoeffs_.set
(
patchI,
new Field<Type>
(
psi.mesh().boundary()[patchI].size(),
Zero
)
);
if (tmat.movable())
{
faceFluxCorrectionPtr_ = tmat().faceFluxCorrectionPtr_;
tmat().faceFluxCorrectionPtr_ = nullptr;
}
else
{
faceFluxCorrectionPtr_ =
new GeometricField<Type, faePatchField, edgeMesh>
(
*(tmat().faceFluxCorrectionPtr_)
);
}
}
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::faMatrix<Type>::clone() const
{
return tmp<faMatrix<Type>>::New(*this);
tmat.clear();
}
// * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * * //
template<class Type>
Foam::faMatrix<Type>::~faMatrix()
{
@ -601,7 +573,7 @@ void Foam::faMatrix<Type>::relax(const scalar alpha)
}
// Finally add the relaxation contribution to the source.
S += (D - D0)*psi_.internalField();
S += (D - D0)*psi_.primitiveField();
}
@ -690,7 +662,7 @@ Foam::faMatrix<Type>::H() const
Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt);
}
Hphi.primitiveFieldRef() += lduMatrix::H(psi_.internalField()) + source_;
Hphi.primitiveFieldRef() += lduMatrix::H(psi_.primitiveField()) + source_;
addBoundarySource(Hphi.primitiveFieldRef());
Hphi.primitiveFieldRef() /= psi_.mesh().S();
@ -908,11 +880,22 @@ void Foam::faMatrix<Type>::operator-=(const tmp<faMatrix<Type>>& tfamv)
template<class Type>
void Foam::faMatrix<Type>::operator+=
(
const GeometricField<Type, faPatchField, areaMesh>& su
const DimensionedField<Type, areaMesh>& su
)
{
checkMethod(*this, su, "+=");
source() -= su.mesh().S()*su.internalField();
source() -= su.mesh().S()*su.field();
}
template<class Type>
void Foam::faMatrix<Type>::operator+=
(
const tmp<DimensionedField<Type, areaMesh>>& tsu
)
{
operator+=(tsu());
tsu.clear();
}
@ -930,11 +913,22 @@ void Foam::faMatrix<Type>::operator+=
template<class Type>
void Foam::faMatrix<Type>::operator-=
(
const GeometricField<Type, faPatchField, areaMesh>& su
const DimensionedField<Type, areaMesh>& su
)
{
checkMethod(*this, su, "-=");
source() += su.mesh().S()*su.internalField();
source() += su.mesh().S()*su.field();
}
template<class Type>
void Foam::faMatrix<Type>::operator-=
(
const tmp<DimensionedField<Type, areaMesh>>& tsu
)
{
operator-=(tsu());
tsu.clear();
}
@ -955,7 +949,7 @@ void Foam::faMatrix<Type>::operator+=
const dimensioned<Type>& su
)
{
source() -= su.mesh().S()*su;
source() -= psi().mesh().S()*su;
}
@ -965,25 +959,38 @@ void Foam::faMatrix<Type>::operator-=
const dimensioned<Type>& su
)
{
source() += su.mesh().S()*su;
source() += psi().mesh().S()*su;
}
template<class Type>
void Foam::faMatrix<Type>::operator+=(const Foam::zero)
{}
template<class Type>
void Foam::faMatrix<Type>::operator-=(const Foam::zero)
{}
template<class Type>
void Foam::faMatrix<Type>::operator*=
(
const areaScalarField& vsf
const areaScalarField::Internal& dsf
)
{
dimensions_ *= vsf.dimensions();
lduMatrix::operator*=(vsf.internalField());
source_ *= vsf.internalField();
dimensions_ *= dsf.dimensions();
lduMatrix::operator*=(dsf.field());
source_ *= dsf.field();
forAll(boundaryCoeffs_, patchI)
forAll(boundaryCoeffs_, patchi)
{
const scalarField psf(vsf.boundaryField()[patchI].patchInternalField());
internalCoeffs_[patchI] *= psf;
boundaryCoeffs_[patchI] *= psf;
const scalarField pisf
(
dsf.mesh().boundary()[patchi].patchInternalField(dsf.field())
);
internalCoeffs_[patchi] *= pisf;
boundaryCoeffs_[patchi] *= pisf;
}
if (faceFluxCorrectionPtr_)
@ -998,11 +1005,22 @@ void Foam::faMatrix<Type>::operator*=
template<class Type>
void Foam::faMatrix<Type>::operator*=
(
const tmp<areaScalarField>& tvsf
const tmp<areaScalarField::Internal>& tfld
)
{
operator*=(tvsf());
tvsf.clear();
operator*=(tfld());
tfld.clear();
}
template<class Type>
void Foam::faMatrix<Type>::operator*=
(
const tmp<areaScalarField>& tfld
)
{
operator*=(tfld());
tfld.clear();
}
@ -1030,34 +1048,32 @@ void Foam::faMatrix<Type>::operator*=
template<class Type>
void Foam::checkMethod
(
const faMatrix<Type>& fam1,
const faMatrix<Type>& fam2,
const faMatrix<Type>& mat1,
const faMatrix<Type>& mat2,
const char* op
)
{
if (&fam1.psi() != &fam2.psi())
if (&mat1.psi() != &mat2.psi())
{
FatalErrorInFunction
<< "incompatible fields for operation "
<< endl << " "
<< "[" << fam1.psi().name() << "] "
<< "Incompatible fields for operation\n "
<< "[" << mat1.psi().name() << "] "
<< op
<< " [" << fam2.psi().name() << "]"
<< " [" << mat2.psi().name() << "]"
<< abort(FatalError);
}
if
(
dimensionSet::checking()
&& fam1.dimensions() != fam2.dimensions()
&& mat1.dimensions() != mat2.dimensions()
)
{
FatalErrorInFunction
<< "incompatible dimensions for operation "
<< endl << " "
<< "[" << fam1.psi().name() << fam1.dimensions()/dimArea << " ] "
<< "Incompatible dimensions for operation\n "
<< "[" << mat1.psi().name() << mat1.dimensions()/dimArea << " ] "
<< op
<< " [" << fam2.psi().name() << fam2.dimensions()/dimArea << " ]"
<< " [" << mat2.psi().name() << mat2.dimensions()/dimArea << " ]"
<< abort(FatalError);
}
}
@ -1066,23 +1082,22 @@ void Foam::checkMethod
template<class Type>
void Foam::checkMethod
(
const faMatrix<Type>& fam,
const GeometricField<Type, faPatchField, areaMesh>& vf,
const faMatrix<Type>& mat,
const DimensionedField<Type, areaMesh>& fld,
const char* op
)
{
if
(
dimensionSet::checking()
&& fam.dimensions()/dimArea != vf.dimensions()
&& mat.dimensions()/dimArea != fld.dimensions()
)
{
FatalErrorInFunction
<< "incompatible dimensions for operation "
<< endl << " "
<< "[" << fam.psi().name() << fam.dimensions()/dimArea << " ] "
<< "Incompatible dimensions for operation\n "
<< "[" << mat.psi().name() << mat.dimensions()/dimArea << " ] "
<< op
<< " [" << vf.name() << vf.dimensions() << " ]"
<< " [" << fld.name() << fld.dimensions() << " ]"
<< abort(FatalError);
}
}
@ -1091,7 +1106,7 @@ void Foam::checkMethod
template<class Type>
void Foam::checkMethod
(
const faMatrix<Type>& fam,
const faMatrix<Type>& mat,
const dimensioned<Type>& dt,
const char* op
)
@ -1099,13 +1114,12 @@ void Foam::checkMethod
if
(
dimensionSet::checking()
&& fam.dimensions()/dimArea != dt.dimensions()
&& mat.dimensions()/dimArea != dt.dimensions()
)
{
FatalErrorInFunction
<< "incompatible dimensions for operation "
<< endl << " "
<< "[" << fam.psi().name() << fam.dimensions()/dimArea << " ] "
<< "Incompatible dimensions for operation\n "
<< "[" << mat.psi().name() << mat.dimensions()/dimArea << " ] "
<< op
<< " [" << dt.name() << dt.dimensions() << " ]"
<< abort(FatalError);
@ -1117,7 +1131,7 @@ template<class Type>
Foam::SolverPerformance<Type> Foam::solve
(
faMatrix<Type>& fam,
Istream& solverControls
const dictionary& solverControls
)
{
return fam.solve(solverControls);
@ -1127,14 +1141,14 @@ Foam::SolverPerformance<Type> Foam::solve
template<class Type>
Foam::SolverPerformance<Type> Foam::solve
(
const tmp<faMatrix<Type>>& tfam,
Istream& solverControls
const tmp<faMatrix<Type>>& tmat,
const dictionary& solverControls
)
{
SolverPerformance<Type> solverPerf =
const_cast<faMatrix<Type>&>(tfam()).solve(solverControls);
SolverPerformance<Type> solverPerf(tmat.constCast().solve(solverControls));
tmat.clear();
tfam.clear();
return solverPerf;
}
@ -1147,12 +1161,12 @@ Foam::SolverPerformance<Type> Foam::solve(faMatrix<Type>& fam)
template<class Type>
Foam::SolverPerformance<Type> Foam::solve(const tmp<faMatrix<Type>>& tfam)
Foam::SolverPerformance<Type> Foam::solve(const tmp<faMatrix<Type>>& tmat)
{
SolverPerformance<Type> solverPerf =
const_cast<faMatrix<Type>&>(tfam()).solve();
SolverPerformance<Type> solverPerf(tmat.constCast().solve());
tmat.clear();
tfam.clear();
return solverPerf;
}
@ -1350,12 +1364,12 @@ template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
(
const faMatrix<Type>& A,
const GeometricField<Type, faPatchField, areaMesh>& su
const DimensionedField<Type, areaMesh>& su
)
{
checkMethod(A, su, "+");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
tC.ref().source() -= su.mesh().S()*su.internalField();
auto tC(A.clone());
tC.ref().source() -= su.mesh().S()*su.field();
return tC;
}
@ -1363,13 +1377,14 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
(
const tmp<faMatrix<Type>>& tA,
const GeometricField<Type, faPatchField, areaMesh>& su
const faMatrix<Type>& A,
const tmp<DimensionedField<Type, areaMesh>>& tsu
)
{
checkMethod(tA(), su, "+");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() -= su.mesh().S()*su.internalField();
checkMethod(A, tsu(), "+");
auto tC(A.clone());
tC.ref().source() -= tsu().mesh().S()*tsu().field();
tsu.clear();
return tC;
}
@ -1382,8 +1397,37 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
)
{
checkMethod(A, tsu(), "+");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
auto tC(A.clone());
tC.ref().source() -= tsu().mesh().S()*tsu().primitiveField();
tsu.clear();
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
(
const tmp<faMatrix<Type>>& tA,
const DimensionedField<Type, areaMesh>& su
)
{
checkMethod(tA(), su, "+");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() -= su.mesh().S()*su.field();
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
(
const tmp<faMatrix<Type>>& tA,
const tmp<DimensionedField<Type, areaMesh>>& tsu
)
{
checkMethod(tA(), tsu(), "+");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() -= tsu().mesh().S()*tsu().field();
tsu.clear();
return tC;
}
@ -1398,7 +1442,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
{
checkMethod(tA(), tsu(), "+");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
tC.ref().source() -= tsu().mesh().S()*tsu().primitiveField();
tsu.clear();
return tC;
}
@ -1407,13 +1451,13 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
(
const GeometricField<Type, faPatchField, areaMesh>& su,
const DimensionedField<Type, areaMesh>& su,
const faMatrix<Type>& A
)
{
checkMethod(A, su, "+");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
tC.ref().source() -= su.mesh().S()*su.internalField();
auto tC(A.clone());
tC.ref().source() -= su.mesh().S()*su.field();
return tC;
}
@ -1421,13 +1465,14 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
(
const GeometricField<Type, faPatchField, areaMesh>& su,
const tmp<faMatrix<Type>>& tA
const tmp<DimensionedField<Type, areaMesh>>& tsu,
const faMatrix<Type>& A
)
{
checkMethod(tA(), su, "+");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() -= su.mesh().S()*su.internalField();
checkMethod(A, tsu(), "+");
auto tC(A.clone());
tC.ref().source() -= tsu().mesh().S()*tsu().field();
tsu.clear();
return tC;
}
@ -1440,8 +1485,37 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
)
{
checkMethod(A, tsu(), "+");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
auto tC(A.clone());
tC.ref().source() -= tsu().mesh().S()*tsu().primitiveField();
tsu.clear();
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
(
const DimensionedField<Type, areaMesh>& su,
const tmp<faMatrix<Type>>& tA
)
{
checkMethod(tA(), su, "+");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() -= su.mesh().S()*su.field();
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
(
const tmp<DimensionedField<Type, areaMesh>>& tsu,
const tmp<faMatrix<Type>>& tA
)
{
checkMethod(tA(), tsu(), "+");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() -= tsu().mesh().S()*tsu().field();
tsu.clear();
return tC;
}
@ -1456,7 +1530,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
{
checkMethod(tA(), tsu(), "+");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
tC.ref().source() -= tsu().mesh().S()*tsu().primitiveField();
tsu.clear();
return tC;
}
@ -1466,12 +1540,12 @@ template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
(
const faMatrix<Type>& A,
const GeometricField<Type, faPatchField, areaMesh>& su
const DimensionedField<Type, areaMesh>& su
)
{
checkMethod(A, su, "-");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
tC.ref().source() += su.mesh().S()*su.internalField();
auto tC(A.clone());
tC.ref().source() += su.mesh().S()*su.field();
return tC;
}
@ -1479,13 +1553,14 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
(
const tmp<faMatrix<Type>>& tA,
const GeometricField<Type, faPatchField, areaMesh>& su
const faMatrix<Type>& A,
const tmp<DimensionedField<Type, areaMesh>>& tsu
)
{
checkMethod(tA(), su, "-");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() += su.mesh().S()*su.internalField();
checkMethod(A, tsu(), "-");
auto tC(A.clone());
tC.ref().source() += tsu().mesh().S()*tsu().field();
tsu.clear();
return tC;
}
@ -1498,8 +1573,37 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
)
{
checkMethod(A, tsu(), "-");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
tC.ref().source() += tsu().mesh().S()*tsu().internalField();
auto tC(A.clone());
tC.ref().source() += tsu().mesh().S()*tsu().primitiveField();
tsu.clear();
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
(
const tmp<faMatrix<Type>>& tA,
const DimensionedField<Type, areaMesh>& su
)
{
checkMethod(tA(), su, "-");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() += su.mesh().S()*su.field();
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
(
const tmp<faMatrix<Type>>& tA,
const tmp<DimensionedField<Type, areaMesh>>& tsu
)
{
checkMethod(tA(), tsu(), "-");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() += tsu().mesh().S()*tsu().field();
tsu.clear();
return tC;
}
@ -1514,7 +1618,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
{
checkMethod(tA(), tsu(), "-");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() += tsu().mesh().S()*tsu().internalField();
tC.ref().source() += tsu().mesh().S()*tsu().primitiveField();
tsu.clear();
return tC;
}
@ -1523,14 +1627,14 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
(
const GeometricField<Type, faPatchField, areaMesh>& su,
const DimensionedField<Type, areaMesh>& su,
const faMatrix<Type>& A
)
{
checkMethod(A, su, "-");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
auto tC(A.clone());
tC.ref().negate();
tC.ref().source() -= su.mesh().S()*su.internalField();
tC.ref().source() -= su.mesh().S()*su.field();
return tC;
}
@ -1538,14 +1642,15 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
(
const GeometricField<Type, faPatchField, areaMesh>& su,
const tmp<faMatrix<Type>>& tA
const tmp<DimensionedField<Type, areaMesh>>& tsu,
const faMatrix<Type>& A
)
{
checkMethod(tA(), su, "-");
tmp<faMatrix<Type>> tC(tA.ptr());
checkMethod(A, tsu(), "-");
auto tC(A.clone());
tC.ref().negate();
tC.ref().source() -= su.mesh().S()*su.internalField();
tC.ref().source() -= tsu().mesh().S()*tsu().field();
tsu.clear();
return tC;
}
@ -1558,9 +1663,40 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
)
{
checkMethod(A, tsu(), "-");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
auto tC(A.clone());
tC.ref().negate();
tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
tC.ref().source() -= tsu().mesh().S()*tsu().primitiveField();
tsu.clear();
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
(
const DimensionedField<Type, areaMesh>& su,
const tmp<faMatrix<Type>>& tA
)
{
checkMethod(tA(), su, "-");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().negate();
tC.ref().source() -= su.mesh().S()*su.field();
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
(
const tmp<DimensionedField<Type, areaMesh>>& tsu,
const tmp<faMatrix<Type>>& tA
)
{
checkMethod(tA(), tsu(), "-");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().negate();
tC.ref().source() -= tsu().mesh().S()*tsu().field();
tsu.clear();
return tC;
}
@ -1576,7 +1712,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
checkMethod(tA(), tsu(), "-");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().negate();
tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
tC.ref().source() -= tsu().mesh().S()*tsu().primitiveField();
tsu.clear();
return tC;
}
@ -1590,7 +1726,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
)
{
checkMethod(A, su, "+");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
auto tC(A.clone());
tC.ref().source() -= su.value()*A.psi().mesh().S();
return tC;
}
@ -1618,7 +1754,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator+
)
{
checkMethod(A, su, "+");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
auto tC(A.clone());
tC.ref().source() -= su.value()*A.psi().mesh().S();
return tC;
}
@ -1646,7 +1782,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
)
{
checkMethod(A, su, "-");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
auto tC(A.clone());
tC.ref().source() += su.value()*tC().psi().mesh().S();
return tC;
}
@ -1674,7 +1810,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator-
)
{
checkMethod(A, su, "-");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
auto tC(A.clone());
tC.ref().negate();
tC.ref().source() -= su.value()*A.psi().mesh().S();
return tC;
@ -1700,12 +1836,12 @@ template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
(
const faMatrix<Type>& A,
const GeometricField<Type, faPatchField, areaMesh>& su
const DimensionedField<Type, areaMesh>& su
)
{
checkMethod(A, su, "==");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
tC.ref().source() += su.mesh().S()*su.internalField();
auto tC(A.clone());
tC.ref().source() += su.mesh().S()*su.field();
return tC;
}
@ -1713,13 +1849,14 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
(
const tmp<faMatrix<Type>>& tA,
const GeometricField<Type, faPatchField, areaMesh>& su
const faMatrix<Type>& A,
const tmp<DimensionedField<Type, areaMesh>>& tsu
)
{
checkMethod(tA(), su, "==");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() += su.mesh().S()*su.internalField();
checkMethod(A, tsu(), "==");
auto tC(A.clone());
tC.ref().source() += tsu().mesh().S()*tsu().field();
tsu.clear();
return tC;
}
@ -1732,8 +1869,37 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
)
{
checkMethod(A, tsu(), "==");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
tC.ref().source() += tsu().mesh().S()*tsu().internalField();
auto tC(A.clone());
tC.ref().source() += tsu().mesh().S()*tsu().primitiveField();
tsu.clear();
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
(
const tmp<faMatrix<Type>>& tA,
const DimensionedField<Type, areaMesh>& su
)
{
checkMethod(tA(), su, "==");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() += su.mesh().S()*su.field();
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
(
const tmp<faMatrix<Type>>& tA,
const tmp<DimensionedField<Type, areaMesh>>& tsu
)
{
checkMethod(tA(), tsu(), "==");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() += tsu().mesh().S()*tsu().field();
tsu.clear();
return tC;
}
@ -1748,7 +1914,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
{
checkMethod(tA(), tsu(), "==");
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref().source() += tsu().mesh().S()*tsu().internalField();
tC.ref().source() += tsu().mesh().S()*tsu().primitiveField();
tsu.clear();
return tC;
}
@ -1762,7 +1928,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
)
{
checkMethod(A, su, "==");
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
auto tC(A.clone());
tC.ref().source() += A.psi().mesh().S()*su.value();
return tC;
}
@ -1782,15 +1948,50 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
(
const faMatrix<Type>& A,
const Foam::zero
)
{
return A;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator==
(
const tmp<faMatrix<Type>>& tA,
const Foam::zero
)
{
return tA;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator*
(
const areaScalarField& vsf,
const areaScalarField::Internal& dsf,
const faMatrix<Type>& A
)
{
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
tC.ref() *= vsf;
auto tC(A.clone());
tC.ref() *= dsf;
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator*
(
const tmp<areaScalarField::Internal>& tdsf,
const faMatrix<Type>& A
)
{
auto tC(A.clone());
tC.ref() *= tdsf;
return tC;
}
@ -1802,7 +2003,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator*
const faMatrix<Type>& A
)
{
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
auto tC(A.clone());
tC.ref() *= tvsf;
return tC;
}
@ -1811,12 +2012,25 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator*
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator*
(
const areaScalarField& vsf,
const areaScalarField::Internal& dsf,
const tmp<faMatrix<Type>>& tA
)
{
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref() *= vsf;
tC.ref() *= dsf;
return tC;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>> Foam::operator*
(
const tmp<areaScalarField::Internal>& tdsf,
const tmp<faMatrix<Type>>& tA
)
{
tmp<faMatrix<Type>> tC(tA.ptr());
tC.ref() *= tdsf;
return tC;
}
@ -1841,7 +2055,7 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator*
const faMatrix<Type>& A
)
{
tmp<faMatrix<Type>> tC(new faMatrix<Type>(A));
auto tC(A.clone());
tC.ref() *= ds;
return tC;
}

225
src/finiteArea/faMatrices/faMatrix/faMatrix.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2020-2021 OpenCFD Ltd.
Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -40,8 +40,8 @@ Author
\*---------------------------------------------------------------------------*/
#ifndef faMatrix_H
#define faMatrix_H
#ifndef Foam_faMatrix_H
#define Foam_faMatrix_H
#include "areaFields.H"
#include "edgeFields.H"
@ -49,6 +49,7 @@ Author
#include "tmp.H"
#include "autoPtr.H"
#include "dimensionedTypes.H"
#include "zero.H"
#include "className.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -78,7 +79,7 @@ public:
// Public Types
//- Field type for psi
//- The geometric field type for psi
typedef
GeometricField<Type, faPatchField, areaMesh>
psiFieldType;
@ -88,6 +89,8 @@ public:
GeometricField<Type, faePatchField, edgeMesh>
faceFluxFieldType;
/// //- The internal field type for the psi field
/// typedef DimensionedField<Type, areaMesh> Internal;
private:
@ -219,6 +222,7 @@ public:
};
// Runtime information
ClassName("faMatrix");
@ -234,15 +238,26 @@ public:
//- Copy construct
faMatrix(const faMatrix<Type>&);
//- Construct from Istream given field to solve for
//- Copy/move construct from tmp<faMatrix<Type>>
faMatrix(const tmp<faMatrix<Type>>&);
//- Deprecated(2022-05) - construct with dimensionSet instead
// \deprecated(2022-05) - construct with dimensionSet instead
FOAM_DEPRECATED_FOR(2022-05, "Construct with dimensionSet")
faMatrix
(
const GeometricField<Type, faPatchField, areaMesh>& psi,
Istream& is
);
)
:
faMatrix<Type>(psi, dimensionSet(is))
{}
//- Clone
tmp<faMatrix<Type>> clone() const;
//- Construct and return a clone
tmp<faMatrix<Type>> clone() const
{
return tmp<faMatrix<Type>>::New(*this);
}
//- Destructor
@ -258,45 +273,45 @@ public:
return psi_;
}
const dimensionSet& dimensions() const
const dimensionSet& dimensions() const noexcept
{
return dimensions_;
}
Field<Type>& source()
Field<Type>& source() noexcept
{
return source_;
}
const Field<Type>& source() const
const Field<Type>& source() const noexcept
{
return source_;
}
//- faBoundary scalar field containing pseudo-matrix coeffs
//- for internal cells
const FieldField<Field, Type>& internalCoeffs() const
const FieldField<Field, Type>& internalCoeffs() const noexcept
{
return internalCoeffs_;
}
//- faBoundary scalar field containing pseudo-matrix coeffs
//- for internal cells
FieldField<Field, Type>& internalCoeffs()
FieldField<Field, Type>& internalCoeffs() noexcept
{
return internalCoeffs_;
}
//- faBoundary scalar field containing pseudo-matrix coeffs
//- for boundary cells
const FieldField<Field, Type>& boundaryCoeffs() const
const FieldField<Field, Type>& boundaryCoeffs() const noexcept
{
return boundaryCoeffs_;
}
//- faBoundary scalar field containing pseudo-matrix coeffs
//- for boundary cells
FieldField<Field, Type>& boundaryCoeffs()
FieldField<Field, Type>& boundaryCoeffs() noexcept
{
return boundaryCoeffs_;
}
@ -390,7 +405,7 @@ public:
void relax();
//- Solve returning the solution statistics.
// Solver controls read Istream
// Use the given solver controls
SolverPerformance<Type> solve(const dictionary&);
//- Solve returning the solution statistics.
@ -418,6 +433,7 @@ public:
void operator=(const faMatrix<Type>&);
void operator=(const tmp<faMatrix<Type>>&);
//- Inplace negate
void negate();
void operator+=(const faMatrix<Type>&);
@ -426,16 +442,28 @@ public:
void operator-=(const faMatrix<Type>&);
void operator-=(const tmp<faMatrix<Type>>&);
void operator+=(const GeometricField<Type,faPatchField,areaMesh>&);
void operator+=(const tmp<GeometricField<Type,faPatchField,areaMesh>>&);
void operator+=(const DimensionedField<Type, areaMesh>&);
void operator+=(const tmp<DimensionedField<Type, areaMesh>>&);
void operator+=
(
const tmp<GeometricField<Type, faPatchField, areaMesh>>&
);
void operator-=(const GeometricField<Type,faPatchField,areaMesh>&);
void operator-=(const tmp<GeometricField<Type,faPatchField,areaMesh>>&);
void operator-=(const DimensionedField<Type, areaMesh>&);
void operator-=(const tmp<DimensionedField<Type, areaMesh>>&);
void operator-=
(
const tmp<GeometricField<Type, faPatchField, areaMesh>>&
);
void operator+=(const dimensioned<Type>&);
void operator-=(const dimensioned<Type>&);
void operator*=(const areaScalarField&);
void operator+=(const Foam::zero);
void operator-=(const Foam::zero);
void operator*=(const areaScalarField::Internal&);
void operator*=(const tmp<areaScalarField::Internal>&);
void operator*=(const tmp<areaScalarField>&);
void operator*=(const dimensioned<scalar>&);
@ -465,7 +493,7 @@ template<class Type>
void checkMethod
(
const faMatrix<Type>&,
const GeometricField<Type, faPatchField, areaMesh>&,
const DimensionedField<Type, areaMesh>&,
const char*
);
@ -479,16 +507,16 @@ void checkMethod
//- Solve returning the solution statistics given convergence tolerance
// Solver controls read Istream
// Use the given solver controls
template<class Type>
SolverPerformance<Type> solve(faMatrix<Type>&, Istream&);
SolverPerformance<Type> solve(faMatrix<Type>&, const dictionary&);
//- Solve returning the solution statistics given convergence tolerance,
//- deleting temporary matrix after solution.
// Solver controls read Istream
// Use the given solver controls
template<class Type>
SolverPerformance<Type> solve(const tmp<faMatrix<Type>>&, Istream&);
SolverPerformance<Type> solve(const tmp<faMatrix<Type>>&, const dictionary&);
//- Solve returning the solution statistics given convergence tolerance
@ -506,18 +534,21 @@ SolverPerformance<Type> solve(const tmp<faMatrix<Type>>&);
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
//- Unary negation
template<class Type>
tmp<faMatrix<Type>> operator-
(
const faMatrix<Type>&
);
//- Unary negation
template<class Type>
tmp<faMatrix<Type>> operator-
(
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator+
(
@ -546,6 +577,7 @@ tmp<faMatrix<Type>> operator+
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator-
(
@ -574,6 +606,7 @@ tmp<faMatrix<Type>> operator-
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator==
(
@ -602,18 +635,19 @@ tmp<faMatrix<Type>> operator==
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator+
(
const faMatrix<Type>&,
const GeometricField<Type, faPatchField, areaMesh>&
const DimensionedField<Type, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> operator+
(
const tmp<faMatrix<Type>>&,
const GeometricField<Type, faPatchField, areaMesh>&
const faMatrix<Type>&,
const tmp<DimensionedField<Type, areaMesh>>&
);
template<class Type>
@ -623,6 +657,20 @@ tmp<faMatrix<Type>> operator+
const tmp<GeometricField<Type, faPatchField, areaMesh>>&
);
template<class Type>
tmp<faMatrix<Type>> operator+
(
const tmp<faMatrix<Type>>&,
const DimensionedField<Type, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> operator+
(
const tmp<faMatrix<Type>>&,
const tmp<DimensionedField<Type, areaMesh>>&
);
template<class Type>
tmp<faMatrix<Type>> operator+
(
@ -633,15 +681,15 @@ tmp<faMatrix<Type>> operator+
template<class Type>
tmp<faMatrix<Type>> operator+
(
const GeometricField<Type, faPatchField, areaMesh>&,
const DimensionedField<Type, areaMesh>&,
const faMatrix<Type>&
);
template<class Type>
tmp<faMatrix<Type>> operator+
(
const GeometricField<Type, faPatchField, areaMesh>&,
const tmp<faMatrix<Type>>&
const tmp<DimensionedField<Type, areaMesh>>&,
const faMatrix<Type>&
);
template<class Type>
@ -651,6 +699,20 @@ tmp<faMatrix<Type>> operator+
const faMatrix<Type>&
);
template<class Type>
tmp<faMatrix<Type>> operator+
(
const DimensionedField<Type, areaMesh>&,
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator+
(
const tmp<DimensionedField<Type, areaMesh>>&,
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator+
(
@ -662,14 +724,14 @@ template<class Type>
tmp<faMatrix<Type>> operator-
(
const faMatrix<Type>&,
const GeometricField<Type, faPatchField, areaMesh>&
const DimensionedField<Type, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> operator-
(
const tmp<faMatrix<Type>>&,
const GeometricField<Type, faPatchField, areaMesh>&
const faMatrix<Type>&,
const tmp<DimensionedField<Type, areaMesh>>&
);
template<class Type>
@ -683,21 +745,36 @@ template<class Type>
tmp<faMatrix<Type>> operator-
(
const tmp<faMatrix<Type>>&,
const tmp<GeometricField<Type, faPatchField, areaMesh>>&
const DimensionedField<Type, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> operator-
(
const GeometricField<Type, faPatchField, areaMesh>&,
const tmp<faMatrix<Type>>&,
const tmp<DimensionedField<Type, areaMesh>>&
);
template<class Type>
tmp<faMatrix<Type>> operator-
(
const tmp<faMatrix<Type>>&,
const tmp<GeometricField<Type, faPatchField, areaMesh>>&
);
template<class Type>
tmp<faMatrix<Type>> operator-
(
const DimensionedField<Type, areaMesh>&,
const faMatrix<Type>&
);
template<class Type>
tmp<faMatrix<Type>> operator-
(
const GeometricField<Type, faPatchField, areaMesh>&,
const tmp<faMatrix<Type>>&
const tmp<DimensionedField<Type, areaMesh>>&,
const faMatrix<Type>&
);
template<class Type>
@ -707,6 +784,20 @@ tmp<faMatrix<Type>> operator-
const faMatrix<Type>&
);
template<class Type>
tmp<faMatrix<Type>> operator-
(
const DimensionedField<Type, areaMesh>&,
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator-
(
const tmp<DimensionedField<Type, areaMesh>>&,
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator-
(
@ -774,14 +865,14 @@ template<class Type>
tmp<faMatrix<Type>> operator==
(
const faMatrix<Type>&,
const GeometricField<Type, faPatchField, areaMesh>&
const DimensionedField<Type, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> operator==
(
const tmp<faMatrix<Type>>&,
const GeometricField<Type, faPatchField, areaMesh>&
const faMatrix<Type>&,
const tmp<DimensionedField<Type, areaMesh>>&
);
template<class Type>
@ -791,6 +882,20 @@ tmp<faMatrix<Type>> operator==
const tmp<GeometricField<Type, faPatchField, areaMesh>>&
);
template<class Type>
tmp<faMatrix<Type>> operator==
(
const tmp<faMatrix<Type>>&,
const DimensionedField<Type, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> operator==
(
const tmp<faMatrix<Type>>&,
const tmp<DimensionedField<Type, areaMesh>>&
);
template<class Type>
tmp<faMatrix<Type>> operator==
(
@ -813,18 +918,33 @@ tmp<faMatrix<Type>> operator==
);
template<class Type>
tmp<faMatrix<Type>> operator==
(
const faMatrix<Type>&,
const Foam::zero
);
template<class Type>
tmp<faMatrix<Type>> operator==
(
const tmp<faMatrix<Type>>&,
const Foam::zero
);
template<class Type>
tmp<faMatrix<Type>> operator*
(
const areaScalarField&,
const areaScalarField::Internal&,
const faMatrix<Type>&
);
template<class Type>
tmp<faMatrix<Type>> operator*
(
const areaScalarField&,
const tmp<faMatrix<Type>>&
const tmp<areaScalarField::Internal>&,
const faMatrix<Type>&
);
template<class Type>
@ -834,6 +954,20 @@ tmp<faMatrix<Type>> operator*
const faMatrix<Type>&
);
template<class Type>
tmp<faMatrix<Type>> operator*
(
const areaScalarField::Internal&,
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator*
(
const tmp<areaScalarField::Internal>&,
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator*
(
@ -841,6 +975,7 @@ tmp<faMatrix<Type>> operator*
const tmp<faMatrix<Type>>&
);
template<class Type>
tmp<faMatrix<Type>> operator*
(

332
src/finiteArea/finiteArea/fam/famSup.C
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -26,183 +27,278 @@ License
\*---------------------------------------------------------------------------*/
#include "areaFields.H"
#include "edgeFields.H"
#include "faMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
tmp<faMatrix<Type>>
Su
Foam::zeroField
Foam::fam::Su
(
const GeometricField<Type, faPatchField, areaMesh>& su,
const GeometricField<Type, faPatchField, areaMesh>& vf
const Foam::zero,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
const faMesh& mesh = vf.mesh();
tmp<faMatrix<Type>> tfam
(
new faMatrix<Type>
(
vf,
dimArea*su.dimensions()
)
);
faMatrix<Type>& fam = tfam.ref();
fam.source() -= mesh.S()*su.internalField();
return tfam;
return zeroField();
}
template<class Type>
tmp<faMatrix<Type>>
Su
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::Su
(
const dimensioned<Type>& su,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
auto tmat = tmp<faMatrix<Type>>::New
(
fld,
dimArea*su.dimensions()
);
auto& mat = tmat.ref();
const auto& domain = fld.mesh().S();
if (magSqr(su.value()) > VSMALL)
{
mat.source() -= domain*su.value();
}
return tmat;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::Su
(
const DimensionedField<Type, areaMesh>& su,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
auto tmat = tmp<faMatrix<Type>>::New
(
fld,
dimArea*su.dimensions()
);
auto& mat = tmat.ref();
const auto& domain = fld.mesh().S();
mat.source() -= domain*su.field();
return tmat;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::Su
(
const tmp<DimensionedField<Type, areaMesh>>& tsu,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
tmp<faMatrix<Type>> tmat = fam::Su(tsu(), fld);
tsu.clear();
return tmat;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::Su
(
const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
const GeometricField<Type, faPatchField, areaMesh>& vf
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
tmp<faMatrix<Type>> tfam = fam::Su(tsu(), vf);
tmp<faMatrix<Type>> tmat = fam::Su(tsu(), fld);
tsu.clear();
return tfam;
return tmat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
tmp<faMatrix<Type>>
Sp
Foam::zeroField
Foam::fam::Sp
(
const areaScalarField& sp,
const GeometricField<Type, faPatchField, areaMesh>& vf
const Foam::zero,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
const faMesh& mesh = vf.mesh();
tmp<faMatrix<Type>> tfam
(
new faMatrix<Type>
(
vf,
dimArea*sp.dimensions()*vf.dimensions()
)
);
faMatrix<Type>& fam = tfam.ref();
fam.diag() += mesh.S()*sp.internalField();
return tfam;
}
template<class Type>
tmp<faMatrix<Type>>
Sp
(
const tmp<areaScalarField>& tsp,
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam = fam::Sp(tsp(), vf);
tsp.clear();
return tfam;
return zeroField();
}
template<class Type>
tmp<faMatrix<Type>>
Sp
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::Sp
(
const dimensionedScalar& sp,
const GeometricField<Type, faPatchField, areaMesh>& vf
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
const faMesh& mesh = vf.mesh();
tmp<faMatrix<Type>> tfam
auto tmat = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
dimArea*sp.dimensions()*vf.dimensions()
)
fld,
dimArea*sp.dimensions()*fld.dimensions()
);
faMatrix<Type>& fam = tfam.ref();
auto& mat = tmat.ref();
const auto& domain = fld.mesh().S();
fam.diag() += mesh.S()*sp.value();
if (mag(sp.value()) > ROOTVSMALL)
{
mat.diag() += domain*sp.value();
}
return tfam;
return tmat;
}
template<class Type>
tmp<faMatrix<Type>>
SuSp
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::Sp
(
const areaScalarField& sp,
const GeometricField<Type, faPatchField, areaMesh>& vf
const DimensionedField<scalar, areaMesh>& sp,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
const faMesh& mesh = vf.mesh();
tmp<faMatrix<Type>> tfam
auto tmat = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
dimArea*sp.dimensions()*vf.dimensions()
)
fld,
dimArea*sp.dimensions()*fld.dimensions()
);
faMatrix<Type>& fam = tfam.ref();
auto& mat = tmat.ref();
const auto& domain = fld.mesh().S();
fam.diag() +=
mesh.S()*max
(
sp.internalField(),
dimensionedScalar("0", sp.dimensions(), Zero)
);
mat.diag() += domain*sp.field();
fam.source() -=
mesh.S()*min
(
sp.internalField(),
dimensionedScalar("0", sp.dimensions(), Zero)
)
*vf.internalField();
return tfam;
return tmat;
}
template<class Type>
tmp<faMatrix<Type>>
SuSp
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::Sp
(
const tmp<DimensionedField<scalar, areaMesh>>& tsp,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
tmp<faMatrix<Type>> tmat = fam::Sp(tsp(), fld);
tsp.clear();
return tmat;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::Sp
(
const tmp<areaScalarField>& tsp,
const GeometricField<Type, faPatchField, areaMesh>& vf
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
tmp<faMatrix<Type>> tfam = fam::SuSp(tsp(), vf);
tmp<faMatrix<Type>> tmat = fam::Sp(tsp(), fld);
tsp.clear();
return tfam;
return tmat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fam
template<class Type>
Foam::zeroField
Foam::fam::SuSp
(
const Foam::zero,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
return zeroField();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
template<class Type>
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::SuSp
(
const dimensionedScalar& susp,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
auto tmat = tmp<faMatrix<Type>>::New
(
fld,
dimArea*susp.dimensions()*fld.dimensions()
);
auto& mat = tmat.ref();
const auto& domain = fld.mesh().S();
if (susp.value() > ROOTVSMALL)
{
mat.diag() += domain*susp.value();
}
else if (susp.value() < -ROOTVSMALL)
{
mat.source() -= domain*susp.value()*fld.primitiveField();
}
return tmat;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::SuSp
(
const DimensionedField<scalar, areaMesh>& susp,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
auto tmat = tmp<faMatrix<Type>>::New
(
fld,
dimArea*susp.dimensions()*fld.dimensions()
);
auto& mat = tmat.ref();
const auto& domain = fld.mesh().S();
mat.diag() += domain*max(susp.field(), scalar(0));
mat.source() -= domain*min(susp.field(), scalar(0))*fld.primitiveField();
return tmat;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::SuSp
(
const tmp<DimensionedField<scalar, areaMesh>>& tsusp,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
tmp<faMatrix<Type>> tmat = fam::SuSp(tsusp(), fld);
tsusp.clear();
return tmat;
}
template<class Type>
Foam::tmp<Foam::faMatrix<Type>>
Foam::fam::SuSp
(
const tmp<areaScalarField>& tsusp,
const GeometricField<Type, faPatchField, areaMesh>& fld
)
{
tmp<faMatrix<Type>> tmat = fam::SuSp(tsusp(), fld);
tsusp.clear();
return tmat;
}
// ************************************************************************* //

83
src/finiteArea/finiteArea/fam/famSup.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,18 +28,19 @@ InNamespace
Foam::fam
Description
Calculate the matrix for implicit and explicit sources.
Calculate the finiteArea matrix for implicit and explicit sources.
SourceFiles
famSup.C
\*---------------------------------------------------------------------------*/
#ifndef famSup_H
#define famSup_H
#ifndef Foam_famSup_H
#define Foam_famSup_H
#include "areaFieldsFwd.H"
#include "faMatrix.H"
#include "zeroField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -46,17 +48,40 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Namespace fam functions Declaration
Namespace fam Functions
\*---------------------------------------------------------------------------*/
namespace fam
{
// Explicit source
//- A no-op source
template<class Type>
zeroField Su
(
const Foam::zero,
const GeometricField<Type, faPatchField, areaMesh>&
);
//- A uniform source (no-op for small values)
template<class Type>
tmp<faMatrix<Type>> Su
(
const GeometricField<Type, faPatchField, areaMesh>&,
const dimensioned<Type>& su,
const GeometricField<Type, faPatchField, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> Su
(
const DimensionedField<Type, areaMesh>& su,
const GeometricField<Type, faPatchField, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> Su
(
const tmp<DimensionedField<Type, areaMesh>>&,
const GeometricField<Type, faPatchField, areaMesh>&
);
@ -70,35 +95,73 @@ namespace fam
// Implicit source
//- A no-op source
template<class Type>
zeroField Sp
(
const Foam::zero,
const GeometricField<Type, faPatchField, areaMesh>&
);
//- A uniform source (no-op for small values)
template<class Type>
tmp<faMatrix<Type>> Sp
(
const areaScalarField&,
const dimensionedScalar& sp,
const GeometricField<Type, faPatchField, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> Sp
(
const tmp<areaScalarField>&,
const DimensionedField<scalar, areaMesh>& sp,
const GeometricField<Type, faPatchField, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> Sp
(
const dimensionedScalar&,
const tmp<DimensionedField<scalar, areaMesh>>&,
const GeometricField<Type, faPatchField, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> Sp
(
const tmp<GeometricField<scalar, faPatchField, areaMesh>>&,
const GeometricField<Type, faPatchField, areaMesh>&
);
// Implicit/Explicit source depending on sign of coefficient
//- A no-op source
template<class Type>
zeroField SuSp
(
const Foam::zero,
const GeometricField<Type, faPatchField, areaMesh>&
);
//- A uniform source (no-op for small values)
template<class Type>
tmp<faMatrix<Type>> SuSp
(
const areaScalarField&,
const dimensionedScalar& susp,
const GeometricField<Type, faPatchField, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> SuSp
(
const DimensionedField<scalar, areaMesh>& susp,
const GeometricField<Type, faPatchField, areaMesh>&
);
template<class Type>
tmp<faMatrix<Type>> SuSp
(
const tmp<DimensionedField<scalar, areaMesh>>&,
const GeometricField<Type, faPatchField, areaMesh>&
);

8
src/finiteVolume/cfdTools/general/fvOptions/fvOptionList.H
View File

@ -37,8 +37,8 @@ SourceFile
\*---------------------------------------------------------------------------*/
#ifndef fvOptionList_H
#define fvOptionList_H
#ifndef Foam_fvOptionList_H
#define Foam_fvOptionList_H
#include "fvOption.H"
#include "PtrList.H"
@ -51,7 +51,7 @@ SourceFile
namespace Foam
{
// Forward declaration of friend functions and operators
// Forward Declarations
namespace fv
{
@ -73,7 +73,7 @@ class optionList
{
protected:
// Protected data
// Protected Data
//- Reference to the mesh database
const fvMesh& mesh_;

343
src/finiteVolume/finiteVolume/fvm/fvmSup.C
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -26,71 +27,15 @@ License
\*---------------------------------------------------------------------------*/
#include "volFields.H"
#include "surfaceFields.H"
#include "fvMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Su
(
const DimensionedField<Type, volMesh>& su,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const fvMesh& mesh = vf.mesh();
tmp<fvMatrix<Type>> tfvm
(
new fvMatrix<Type>
(
vf,
dimVol*su.dimensions()
)
);
fvMatrix<Type>& fvm = tfvm.ref();
fvm.source() -= mesh.V()*su.field();
return tfvm;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Su
(
const tmp<DimensionedField<Type, volMesh>>& tsu,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
tmp<fvMatrix<Type>> tfvm = fvm::Su(tsu(), vf);
tsu.clear();
return tfvm;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Su
(
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
tmp<fvMatrix<Type>> tfvm = fvm::Su(tsu(), vf);
tsu.clear();
return tfvm;
}
template<class Type>
Foam::zeroField
Foam::fvm::Su
(
const zero&,
const GeometricField<Type, fvPatchField, volMesh>& vf
const Foam::zero,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
return zeroField();
@ -99,89 +44,86 @@ Foam::fvm::Su
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Sp
Foam::fvm::Su
(
const volScalarField::Internal& sp,
const GeometricField<Type, fvPatchField, volMesh>& vf
const dimensioned<Type>& su,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
const fvMesh& mesh = vf.mesh();
tmp<fvMatrix<Type>> tfvm
auto tmat = tmp<fvMatrix<Type>>::New
(
new fvMatrix<Type>
(
vf,
dimVol*sp.dimensions()*vf.dimensions()
)
fld,
dimVol*su.dimensions()
);
fvMatrix<Type>& fvm = tfvm.ref();
auto& mat = tmat.ref();
const auto& domain = fld.mesh().V();
fvm.diag() += mesh.V()*sp.field();
if (magSqr(su.value()) > VSMALL)
{
mat.source() -= domain*su.value();
}
return tfvm;
return tmat;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Sp
Foam::fvm::Su
(
const tmp<volScalarField::Internal>& tsp,
const GeometricField<Type, fvPatchField, volMesh>& vf
const DimensionedField<Type, volMesh>& su,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
tmp<fvMatrix<Type>> tfvm = fvm::Sp(tsp(), vf);
tsp.clear();
return tfvm;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Sp
(
const tmp<volScalarField>& tsp,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
tmp<fvMatrix<Type>> tfvm = fvm::Sp(tsp(), vf);
tsp.clear();
return tfvm;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Sp
(
const dimensionedScalar& sp,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
const fvMesh& mesh = vf.mesh();
tmp<fvMatrix<Type>> tfvm
auto tmat = tmp<fvMatrix<Type>>::New
(
new fvMatrix<Type>
(
vf,
dimVol*sp.dimensions()*vf.dimensions()
)
fld,
dimVol*su.dimensions()
);
fvMatrix<Type>& fvm = tfvm.ref();
auto& mat = tmat.ref();
const auto& domain = fld.mesh().V();
fvm.diag() += mesh.V()*sp.value();
mat.source() -= domain*su.field();
return tfvm;
return tmat;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Su
(
const tmp<DimensionedField<Type, volMesh>>& tsu,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
tmp<fvMatrix<Type>> tmat = fvm::Su(tsu(), fld);
tsu.clear();
return tmat;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Su
(
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
tmp<fvMatrix<Type>> tmat = fvm::Su(tsu(), fld);
tsu.clear();
return tmat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
Foam::zeroField
Foam::fvm::Sp
(
const zero&,
const Foam::zero,
const GeometricField<Type, fvPatchField, volMesh>&
)
{
@ -191,30 +133,90 @@ Foam::fvm::Sp
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::SuSp
Foam::fvm::Sp
(
const volScalarField::Internal& susp,
const GeometricField<Type, fvPatchField, volMesh>& vf
const dimensionedScalar& sp,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
const fvMesh& mesh = vf.mesh();
tmp<fvMatrix<Type>> tfvm
auto tmat = tmp<fvMatrix<Type>>::New
(
new fvMatrix<Type>
(
vf,
dimVol*susp.dimensions()*vf.dimensions()
)
fld,
dimVol*sp.dimensions()*fld.dimensions()
);
fvMatrix<Type>& fvm = tfvm.ref();
auto& mat = tmat.ref();
const auto& domain = fld.mesh().V();
fvm.diag() += mesh.V()*max(susp.field(), scalar(0));
if (mag(sp.value()) > ROOTVSMALL)
{
mat.diag() += domain*sp.value();
}
fvm.source() -= mesh.V()*min(susp.field(), scalar(0))
*vf.primitiveField();
return tmat;
}
return tfvm;
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Sp
(
const DimensionedField<scalar, volMesh>& sp,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
auto tmat = tmp<fvMatrix<Type>>::New
(
fld,
dimVol*sp.dimensions()*fld.dimensions()
);
auto& mat = tmat.ref();
const auto& domain = fld.mesh().V();
mat.diag() += domain*sp.field();
return tmat;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Sp
(
const tmp<DimensionedField<scalar, volMesh>>& tsp,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
tmp<fvMatrix<Type>> tmat = fvm::Sp(tsp(), fld);
tsp.clear();
return tmat;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Sp
(
const tmp<volScalarField>& tsp,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
tmp<fvMatrix<Type>> tmat = fvm::Sp(tsp(), fld);
tsp.clear();
return tmat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
Foam::zeroField
Foam::fvm::SuSp
(
const Foam::zero,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
return zeroField();
}
@ -222,13 +224,66 @@ template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::SuSp
(
const tmp<volScalarField::Internal>& tsusp,
const GeometricField<Type, fvPatchField, volMesh>& vf
const dimensionedScalar& susp,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
tmp<fvMatrix<Type>> tfvm = fvm::SuSp(tsusp(), vf);
auto tmat = tmp<fvMatrix<Type>>::New
(
fld,
dimVol*susp.dimensions()*fld.dimensions()
);
auto& mat = tmat.ref();
const auto& domain = fld.mesh().V();
if (susp.value() > ROOTVSMALL)
{
mat.diag() += domain*susp.value();
}
else if (susp.value() < -ROOTVSMALL)
{
mat.source() -= domain*susp.value()*fld.primitiveField();
}
return tmat;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::SuSp
(
const DimensionedField<scalar, volMesh>& susp,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
auto tmat = tmp<fvMatrix<Type>>::New
(
fld,
dimVol*susp.dimensions()*fld.dimensions()
);
auto& mat = tmat.ref();
const auto& domain = fld.mesh().V();
mat.diag() += domain*max(susp.field(), scalar(0));
mat.source() -= domain*min(susp.field(), scalar(0))*fld.primitiveField();
return tmat;
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::SuSp
(
const tmp<DimensionedField<scalar, volMesh>>& tsusp,
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
tmp<fvMatrix<Type>> tmat = fvm::SuSp(tsusp(), fld);
tsusp.clear();
return tfvm;
return tmat;
}
@ -237,24 +292,12 @@ Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::SuSp
(
const tmp<volScalarField>& tsusp,
const GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
tmp<fvMatrix<Type>> tfvm = fvm::SuSp(tsusp(), vf);
tmp<fvMatrix<Type>> tmat = fvm::SuSp(tsusp(), fld);
tsusp.clear();
return tfvm;
}
template<class Type>
Foam::zeroField
Foam::fvm::SuSp
(
const zero&,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return zeroField();
return tmat;
}

97
src/finiteVolume/finiteVolume/fvm/fvmSup.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,15 +28,15 @@ InNamespace
Foam::fvm
Description
Calculate the matrix for implicit and explicit sources.
Calculate the finiteVolume matrix for implicit and explicit sources.
SourceFiles
fvmSup.C
\*---------------------------------------------------------------------------*/
#ifndef fvmSup_H
#define fvmSup_H
#ifndef Foam_fvmSup_H
#define Foam_fvmSup_H
#include "volFieldsFwd.H"
#include "fvMatrix.H"
@ -47,17 +48,33 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Namespace fvm functions Declaration
Namespace fvm Functions
\*---------------------------------------------------------------------------*/
namespace fvm
{
// Explicit source
//- A no-op source
template<class Type>
zeroField Su
(
const Foam::zero,
const GeometricField<Type, fvPatchField, volMesh>&
);
//- A uniform source (no-op for small values)
template<class Type>
tmp<fvMatrix<Type>> Su
(
const DimensionedField<Type, volMesh>&,
const dimensioned<Type>& su,
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
tmp<fvMatrix<Type>> Su
(
const DimensionedField<Type, volMesh>& su,
const GeometricField<Type, fvPatchField, volMesh>&
);
@ -75,27 +92,36 @@ namespace fvm
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
zeroField Su
(
const zero&,
const GeometricField<Type, fvPatchField, volMesh>&
);
// Implicit source
//- A no-op source
template<class Type>
zeroField Sp
(
const Foam::zero,
const GeometricField<Type, fvPatchField, volMesh>&
);
//- A uniform source (no-op for small values)
template<class Type>
tmp<fvMatrix<Type>> Sp
(
const volScalarField::Internal&,
const dimensionedScalar& sp,
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
tmp<fvMatrix<Type>> Sp
(
const tmp<volScalarField::Internal>&,
const DimensionedField<scalar, volMesh>& sp,
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
tmp<fvMatrix<Type>> Sp
(
const tmp<DimensionedField<scalar, volMesh>>&,
const GeometricField<Type, fvPatchField, volMesh>&
);
@ -107,35 +133,35 @@ namespace fvm
);
template<class Type>
tmp<fvMatrix<Type>> Sp
(
const dimensionedScalar&,
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
zeroField Sp
(
const zero&,
const GeometricField<Type, fvPatchField, volMesh>&
);
// Implicit/Explicit source depending on sign of coefficient
//- A no-op source
template<class Type>
zeroField SuSp
(
const Foam::zero,
const GeometricField<Type, fvPatchField, volMesh>&
);
//- A uniform source (no-op for small values)
template<class Type>
tmp<fvMatrix<Type>> SuSp
(
const volScalarField::Internal&,
const dimensionedScalar& susp,
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
tmp<fvMatrix<Type>> SuSp
(
const tmp<volScalarField::Internal>&,
const DimensionedField<scalar, volMesh>& susp,
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
tmp<fvMatrix<Type>> SuSp
(
const tmp<DimensionedField<scalar, volMesh>>&,
const GeometricField<Type, fvPatchField, volMesh>&
);
@ -145,13 +171,6 @@ namespace fvm
const tmp<volScalarField>&,
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
zeroField SuSp
(
const zero&,
const GeometricField<Type, fvPatchField, volMesh>&
);
}

7
src/finiteVolume/fvMatrices/fvMatricesFwd.H
View File

@ -28,8 +28,8 @@ Description
\*---------------------------------------------------------------------------*/
#ifndef fvMatricesFwd_H
#define fvMatricesFwd_H
#ifndef Foam_fvMatricesFwd_H
#define Foam_fvMatricesFwd_H
#include "fieldTypes.H"
@ -40,8 +40,7 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
class fvMatrix;
template<class Type> class fvMatrix;
typedef fvMatrix<scalar> fvScalarMatrix;
typedef fvMatrix<vector> fvVectorMatrix;

222
src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2016-2021 OpenCFD Ltd.
Copyright (C) 2016-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -339,7 +339,6 @@ bool Foam::fvMatrix<Type>::checkImplicit(const label fieldi)
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
Foam::fvMatrix<Type>::fvMatrix
(
@ -350,7 +349,7 @@ Foam::fvMatrix<Type>::fvMatrix
lduMatrix(psi.mesh()),
psi_(psi),
useImplicit_(false),
lduAssemblyName_(word::null),
lduAssemblyName_(),
nMatrix_(0),
dimensions_(ds),
source_(psi.size(), Zero),
@ -368,26 +367,18 @@ Foam::fvMatrix<Type>::fvMatrix
internalCoeffs_.set
(
patchi,
new Field<Type>
(
psi.mesh().boundary()[patchi].size(),
Zero
)
new Field<Type>(psi.mesh().boundary()[patchi].size(), Zero)
);
boundaryCoeffs_.set
(
patchi,
new Field<Type>
(
psi.mesh().boundary()[patchi].size(),
Zero
)
new Field<Type>(psi.mesh().boundary()[patchi].size(), Zero)
);
}
auto& psiRef = this->psi(0);
label currentStatePsi = psiRef.eventNo();
const label currentStatePsi = psiRef.eventNo();
psiRef.boundaryFieldRef().updateCoeffs();
psiRef.eventNo() = currentStatePsi;
}
@ -396,7 +387,6 @@ Foam::fvMatrix<Type>::fvMatrix
template<class Type>
Foam::fvMatrix<Type>::fvMatrix(const fvMatrix<Type>& fvm)
:
refCount(),
lduMatrix(fvm),
psi_(fvm.psi_),
useImplicit_(fvm.useImplicit_),
@ -423,113 +413,40 @@ Foam::fvMatrix<Type>::fvMatrix(const fvMatrix<Type>& fvm)
template<class Type>
Foam::fvMatrix<Type>::fvMatrix(const tmp<fvMatrix<Type>>& tfvm)
Foam::fvMatrix<Type>::fvMatrix(const tmp<fvMatrix<Type>>& tmat)
:
lduMatrix
(
const_cast<fvMatrix<Type>&>(tfvm()),
tfvm.isTmp()
),
psi_(tfvm().psi_),
useImplicit_(tfvm().useImplicit_),
lduAssemblyName_(tfvm().lduAssemblyName_),
nMatrix_(tfvm().nMatrix_),
dimensions_(tfvm().dimensions_),
source_
(
const_cast<fvMatrix<Type>&>(tfvm()).source_,
tfvm.isTmp()
),
internalCoeffs_
(
const_cast<fvMatrix<Type>&>(tfvm()).internalCoeffs_,
tfvm.isTmp()
),
boundaryCoeffs_
(
const_cast<fvMatrix<Type>&>(tfvm()).boundaryCoeffs_,
tfvm.isTmp()
),
lduMatrix(tmat.constCast(), tmat.movable()),
psi_(tmat().psi_),
useImplicit_(tmat().useImplicit_),
lduAssemblyName_(tmat().lduAssemblyName_),
nMatrix_(tmat().nMatrix_),
dimensions_(tmat().dimensions_),
source_(tmat.constCast().source_, tmat.movable()),
internalCoeffs_(tmat.constCast().internalCoeffs_, tmat.movable()),
boundaryCoeffs_(tmat.constCast().boundaryCoeffs_, tmat.movable()),
faceFluxCorrectionPtr_(nullptr)
{
DebugInFunction
<< "Copying fvMatrix<Type> for field " << psi_.name() << endl;
<< "Copy/move fvMatrix<Type> for field " << psi_.name() << endl;
if (tfvm().faceFluxCorrectionPtr_)
if (tmat().faceFluxCorrectionPtr_)
{
if (tfvm.isTmp())
if (tmat.movable())
{
faceFluxCorrectionPtr_ = tfvm().faceFluxCorrectionPtr_;
tfvm().faceFluxCorrectionPtr_ = nullptr;
faceFluxCorrectionPtr_ = tmat().faceFluxCorrectionPtr_;
tmat().faceFluxCorrectionPtr_ = nullptr;
}
else
{
faceFluxCorrectionPtr_ =
new GeometricField<Type, fvsPatchField, surfaceMesh>
(
*(tfvm().faceFluxCorrectionPtr_)
*(tmat().faceFluxCorrectionPtr_)
);
}
}
tfvm.clear();
}
template<class Type>
Foam::fvMatrix<Type>::fvMatrix
(
const GeometricField<Type, fvPatchField, volMesh>& psi,
Istream& is
)
:
lduMatrix(psi.mesh()),
psi_(psi),
useImplicit_(false),
lduAssemblyName_(word::null),
nMatrix_(0),
dimensions_(is),
source_(is),
internalCoeffs_(psi.mesh().boundary().size()),
boundaryCoeffs_(psi.mesh().boundary().size()),
faceFluxCorrectionPtr_(nullptr)
{
DebugInFunction
<< "Constructing fvMatrix<Type> for field " << psi_.name() << endl;
checkImplicit();
// Initialise coupling coefficients
forAll(psi.mesh().boundary(), patchi)
{
internalCoeffs_.set
(
patchi,
new Field<Type>
(
psi.mesh().boundary()[patchi].size(),
Zero
)
);
boundaryCoeffs_.set
(
patchi,
new Field<Type>
(
psi.mesh().boundary()[patchi].size(),
Zero
)
);
}
}
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>> Foam::fvMatrix<Type>::clone() const
{
return tmp<fvMatrix<Type>>::New(*this);
tmat.clear();
}
@ -769,21 +686,13 @@ void Foam::fvMatrix<Type>::setBounAndInterCoeffs()
internalCoeffs_.set
(
interfaceI,
new Field<Type>
(
psi.mesh().boundary()[patchi].size(),
Zero
)
new Field<Type>(psi.mesh().boundary()[patchi].size(), Zero)
);
boundaryCoeffs_.set
(
interfaceI,
new Field<Type>
(
psi.mesh().boundary()[patchi].size(),
Zero
)
new Field<Type>(psi.mesh().boundary()[patchi].size(), Zero)
);
interfaceI++;
}
@ -1882,18 +1791,12 @@ void Foam::fvMatrix<Type>::operator-=
template<class Type>
void Foam::fvMatrix<Type>::operator+=
(
const zero&
)
void Foam::fvMatrix<Type>::operator+=(const Foam::zero)
{}
template<class Type>
void Foam::fvMatrix<Type>::operator-=
(
const zero&
)
void Foam::fvMatrix<Type>::operator-=(const Foam::zero)
{}
@ -1930,22 +1833,22 @@ void Foam::fvMatrix<Type>::operator*=
template<class Type>
void Foam::fvMatrix<Type>::operator*=
(
const tmp<volScalarField::Internal>& tdsf
const tmp<volScalarField::Internal>& tfld
)
{
operator*=(tdsf());
tdsf.clear();
operator*=(tfld());
tfld.clear();
}
template<class Type>
void Foam::fvMatrix<Type>::operator*=
(
const tmp<volScalarField>& tvsf
const tmp<volScalarField>& tfld
)
{
operator*=(tvsf());
tvsf.clear();
operator*=(tfld());
tfld.clear();
}
@ -1973,34 +1876,32 @@ void Foam::fvMatrix<Type>::operator*=
template<class Type>
void Foam::checkMethod
(
const fvMatrix<Type>& fvm1,
const fvMatrix<Type>& fvm2,
const fvMatrix<Type>& mat1,
const fvMatrix<Type>& mat2,
const char* op
)
{
if (&fvm1.psi() != &fvm2.psi())
if (&mat1.psi() != &mat2.psi())
{
FatalErrorInFunction
<< "incompatible fields for operation "
<< endl << " "
<< "[" << fvm1.psi().name() << "] "
<< "Incompatible fields for operation\n "
<< "[" << mat1.psi().name() << "] "
<< op
<< " [" << fvm2.psi().name() << "]"
<< " [" << mat2.psi().name() << "]"
<< abort(FatalError);
}
if
(
dimensionSet::checking()
&& fvm1.dimensions() != fvm2.dimensions()
&& mat1.dimensions() != mat2.dimensions()
)
{
FatalErrorInFunction
<< "incompatible dimensions for operation "
<< endl << " "
<< "[" << fvm1.psi().name() << fvm1.dimensions()/dimVolume << " ] "
<< "Incompatible dimensions for operation\n "
<< "[" << mat1.psi().name() << mat1.dimensions()/dimVolume << " ] "
<< op
<< " [" << fvm2.psi().name() << fvm2.dimensions()/dimVolume << " ]"
<< " [" << mat2.psi().name() << mat2.dimensions()/dimVolume << " ]"
<< abort(FatalError);
}
}
@ -2009,22 +1910,22 @@ void Foam::checkMethod
template<class Type>
void Foam::checkMethod
(
const fvMatrix<Type>& fvm,
const DimensionedField<Type, volMesh>& df,
const fvMatrix<Type>& mat,
const DimensionedField<Type, volMesh>& fld,
const char* op
)
{
if
(
dimensionSet::checking()
&& fvm.dimensions()/dimVolume != df.dimensions()
&& mat.dimensions()/dimVolume != fld.dimensions()
)
{
FatalErrorInFunction
<< endl << " "
<< "[" << fvm.psi().name() << fvm.dimensions()/dimVolume << " ] "
<< "Incompatible dimensions for operation\n "
<< "[" << mat.psi().name() << mat.dimensions()/dimVolume << " ] "
<< op
<< " [" << df.name() << df.dimensions() << " ]"
<< " [" << fld.name() << fld.dimensions() << " ]"
<< abort(FatalError);
}
}
@ -2033,7 +1934,7 @@ void Foam::checkMethod
template<class Type>
void Foam::checkMethod
(
const fvMatrix<Type>& fvm,
const fvMatrix<Type>& mat,
const dimensioned<Type>& dt,
const char* op
)
@ -2041,13 +1942,12 @@ void Foam::checkMethod
if
(
dimensionSet::checking()
&& fvm.dimensions()/dimVolume != dt.dimensions()
&& mat.dimensions()/dimVolume != dt.dimensions()
)
{
FatalErrorInFunction
<< "incompatible dimensions for operation "
<< endl << " "
<< "[" << fvm.psi().name() << fvm.dimensions()/dimVolume << " ] "
<< "Incompatible dimensions for operation\n "
<< "[" << mat.psi().name() << mat.dimensions()/dimVolume << " ] "
<< op
<< " [" << dt.name() << dt.dimensions() << " ]"
<< abort(FatalError);
@ -2068,14 +1968,13 @@ Foam::SolverPerformance<Type> Foam::solve
template<class Type>
Foam::SolverPerformance<Type> Foam::solve
(
const tmp<fvMatrix<Type>>& tfvm,
const tmp<fvMatrix<Type>>& tmat,
const dictionary& solverControls
)
{
SolverPerformance<Type> solverPerf =
const_cast<fvMatrix<Type>&>(tfvm()).solve(solverControls);
SolverPerformance<Type> solverPerf(tmat.constCast().solve(solverControls));
tfvm.clear();
tmat.clear();
return solverPerf;
}
@ -2088,12 +1987,11 @@ Foam::SolverPerformance<Type> Foam::solve(fvMatrix<Type>& fvm)
}
template<class Type>
Foam::SolverPerformance<Type> Foam::solve(const tmp<fvMatrix<Type>>& tfvm)
Foam::SolverPerformance<Type> Foam::solve(const tmp<fvMatrix<Type>>& tmat)
{
SolverPerformance<Type> solverPerf =
const_cast<fvMatrix<Type>&>(tfvm()).solve();
SolverPerformance<Type> solverPerf(tmat.constCast().solve());
tfvm.clear();
tmat.clear();
return solverPerf;
}
@ -2289,7 +2187,7 @@ template<class Type>
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
(
const fvMatrix<Type>& A,
const zero&
const Foam::zero
)
{
return A;
@ -2300,7 +2198,7 @@ template<class Type>
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
(
const tmp<fvMatrix<Type>>& tA,
const zero&
const Foam::zero
)
{
return tA;

75
src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2016-2021 OpenCFD Ltd.
Copyright (C) 2016-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -40,8 +40,8 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef fvMatrix_H
#define fvMatrix_H
#ifndef Foam_fvMatrix_H
#define Foam_fvMatrix_H
#include "volFields.H"
#include "surfaceFields.H"
@ -110,7 +110,6 @@ tmp<GeometricField<Type, fvPatchField, volMesh>> operator&
);
/*---------------------------------------------------------------------------*\
Class fvMatrix Declaration
\*---------------------------------------------------------------------------*/
@ -125,7 +124,7 @@ public:
// Public Types
//- Field type for psi
//- The geometric field type for psi
typedef
GeometricField<Type, fvPatchField, volMesh>
psiFieldType;
@ -135,6 +134,10 @@ public:
GeometricField<Type, fvsPatchField, surfaceMesh>
faceFluxFieldType;
/// //- The internal field type for the psi field
/// typedef DimensionedField<Type, volMesh> Internal;
private:
// Private Data
@ -286,6 +289,7 @@ public:
};
// Runtime information
ClassName("fvMatrix");
@ -304,15 +308,23 @@ public:
//- Copy/move construct from tmp<fvMatrix<Type>>
fvMatrix(const tmp<fvMatrix<Type>>&);
//- Construct from Istream given field to solve for
//- Deprecated(2022-05) - construct with dimensionSet instead
// \deprecated(2022-05) - construct with dimensionSet instead
FOAM_DEPRECATED_FOR(2022-05, "Construct with dimensionSet")
fvMatrix
(
const GeometricField<Type, fvPatchField, volMesh>& psi,
Istream& is
);
)
:
fvMatrix<Type>(psi, dimensionSet(is))
{}
//- Clone
tmp<fvMatrix<Type>> clone() const;
//- Construct and return a clone
tmp<fvMatrix<Type>> clone() const
{
return tmp<fvMatrix<Type>>::New(*this);
}
//- Destructor
@ -439,45 +451,45 @@ public:
}
const dimensionSet& dimensions() const
const dimensionSet& dimensions() const noexcept
{
return dimensions_;
}
Field<Type>& source()
Field<Type>& source() noexcept
{
return source_;
}
const Field<Type>& source() const
const Field<Type>& source() const noexcept
{
return source_;
}
//- fvBoundary scalar field containing pseudo-matrix coeffs
//- for internal cells
const FieldField<Field, Type>& internalCoeffs() const
const FieldField<Field, Type>& internalCoeffs() const noexcept
{
return internalCoeffs_;
}
//- fvBoundary scalar field containing pseudo-matrix coeffs
//- for internal cells
FieldField<Field, Type>& internalCoeffs()
FieldField<Field, Type>& internalCoeffs() noexcept
{
return internalCoeffs_;
}
//- fvBoundary scalar field containing pseudo-matrix coeffs
//- for boundary cells
const FieldField<Field, Type>& boundaryCoeffs() const
const FieldField<Field, Type>& boundaryCoeffs() const noexcept
{
return boundaryCoeffs_;
}
//- fvBoundary scalar field containing pseudo-matrix coeffs
//- for boundary cells
FieldField<Field, Type>& boundaryCoeffs()
FieldField<Field, Type>& boundaryCoeffs() noexcept
{
return boundaryCoeffs_;
}
@ -639,6 +651,7 @@ public:
void operator=(const fvMatrix<Type>&);
void operator=(const tmp<fvMatrix<Type>>&);
//- Inplace negate
void negate();
void operator+=(const fvMatrix<Type>&);
@ -647,27 +660,15 @@ public:
void operator-=(const fvMatrix<Type>&);
void operator-=(const tmp<fvMatrix<Type>>&);
void operator+=
(
const DimensionedField<Type, volMesh>&
);
void operator+=
(
const tmp<DimensionedField<Type, volMesh>>&
);
void operator+=(const DimensionedField<Type, volMesh>&);
void operator+=(const tmp<DimensionedField<Type, volMesh>>&);
void operator+=
(
const tmp<GeometricField<Type, fvPatchField, volMesh>>&
);
void operator-=
(
const DimensionedField<Type, volMesh>&
);
void operator-=
(
const tmp<DimensionedField<Type, volMesh>>&
);
void operator-=(const DimensionedField<Type, volMesh>&);
void operator-=(const tmp<DimensionedField<Type, volMesh>>&);
void operator-=
(
const tmp<GeometricField<Type, fvPatchField, volMesh>>&
@ -676,8 +677,8 @@ public:
void operator+=(const dimensioned<Type>&);
void operator-=(const dimensioned<Type>&);
void operator+=(const zero&);
void operator-=(const zero&);
void operator+=(const Foam::zero);
void operator-=(const Foam::zero);
void operator*=(const volScalarField::Internal&);
void operator*=(const tmp<volScalarField::Internal>&);
@ -888,23 +889,25 @@ template<class Type>
tmp<fvMatrix<Type>> operator==
(
const fvMatrix<Type>&,
const zero&
const Foam::zero
);
template<class Type>
tmp<fvMatrix<Type>> operator==
(
const tmp<fvMatrix<Type>>&,
const zero&
const Foam::zero
);
//- Unary negation
template<class Type>
tmp<fvMatrix<Type>> operator-
(
const fvMatrix<Type>&
);
//- Unary negation
template<class Type>
tmp<fvMatrix<Type>> operator-
(

2
src/fvOptions/Make/files
View File

@ -5,7 +5,7 @@ interRegionOption/interRegionOption.C
generalSources=sources/general
$(generalSources)/codedSource/codedFvSources.C
$(generalSources)/semiImplicitSource/semiImplicitSource.C
$(generalSources)/semiImplicitSource/semiImplicitSources.C
derivedSources=sources/derived
$(derivedSources)/acousticDampingSource/acousticDampingSource.C

177
src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020-2021 OpenCFD Ltd.
Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -31,6 +31,7 @@ License
#include "fvMatrices.H"
#include "fvmSup.H"
#include "Constant.H"
#include "Tuple2.H"
// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
@ -46,25 +47,35 @@ Foam::fv::SemiImplicitSource<Type>::volumeModeTypeNames_
});
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
void Foam::fv::SemiImplicitSource<Type>::setFieldData(const dictionary& dict)
void Foam::fv::SemiImplicitSource<Type>::setFieldInjectionRates
(
const dictionary& dict
)
{
label count = dict.size();
fieldNames_.resize(count);
Su_.resize(fieldNames_.size());
Sp_.resize(fieldNames_.size());
fieldNames_.resize_nocopy(count);
Su_.resize(count);
Sp_.resize(count);
fv::option::resetApplied();
count = 0;
for (const entry& dEntry : dict)
{
fieldNames_[count] = dEntry.keyword();
const word& fieldName = dEntry.keyword();
if (!dEntry.isDict())
if (dEntry.isDict())
{
const dictionary& subdict = dEntry.dict();
Su_.set(count, Function1<Type>::New("Su", subdict, &mesh_));
Sp_.set(count, Function1<scalar>::New("Sp", subdict, &mesh_));
}
else
{
Tuple2<Type, scalar> injectionRate;
dEntry.readEntry(injectionRate);
@ -88,21 +99,10 @@ void Foam::fv::SemiImplicitSource<Type>::setFieldData(const dictionary& dict)
)
);
}
else
{
const dictionary& Sdict = dEntry.dict();
Su_.set(count, Function1<Type>::New("Su", Sdict, &mesh_));
Sp_.set(count, Function1<scalar>::New("Sp", Sdict, &mesh_));
}
fieldNames_[count] = fieldName;
++count;
}
// Set volume normalisation
if (volumeMode_ == vmAbsolute)
{
VDash_ = V_;
}
}
@ -119,7 +119,7 @@ Foam::fv::SemiImplicitSource<Type>::SemiImplicitSource
:
fv::cellSetOption(name, modelType, dict, mesh),
volumeMode_(vmAbsolute),
VDash_(1.0)
VDash_(1)
{
read(dict);
}
@ -134,51 +134,7 @@ void Foam::fv::SemiImplicitSource<Type>::addSup
const label fieldi
)
{
if (debug)
{
Info<< "SemiImplicitSource<" << pTraits<Type>::typeName
<< ">::addSup for source " << name_ << endl;
}
const GeometricField<Type, fvPatchField, volMesh>& psi = eqn.psi();
typename GeometricField<Type, fvPatchField, volMesh>::Internal Su
(
IOobject
(
name_ + fieldNames_[fieldi] + "Su",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensioned<Type>(eqn.dimensions()/dimVolume, Zero),
false
);
const scalar tmVal = mesh_.time().timeOutputValue();
UIndirectList<Type>(Su, cells_) = Su_[fieldi].value(tmVal)/VDash_;
volScalarField::Internal Sp
(
IOobject
(
name_ + fieldNames_[fieldi] + "Sp",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensioned<scalar>(Su.dimensions()/psi.dimensions(), Zero),
false
);
UIndirectList<scalar>(Sp, cells_) = Sp_[fieldi].value(tmVal)/VDash_;
eqn += Su + fvm::SuSp(Sp, psi);
return this->addSup(volScalarField::null(), eqn, fieldi);
}
@ -196,17 +152,102 @@ void Foam::fv::SemiImplicitSource<Type>::addSup
<< ">::addSup for source " << name_ << endl;
}
return this->addSup(eqn, fieldi);
const GeometricField<Type, fvPatchField, volMesh>& psi = eqn.psi();
const word& fieldName = fieldNames_[fieldi];
const scalar tmVal = mesh_.time().timeOutputValue();
const dimensionSet SuDims(eqn.dimensions()/dimVolume);
const dimensionSet SpDims(SuDims/psi.dimensions());
// Explicit source
{
const dimensioned<Type> SuValue
(
"Su",
SuDims,
Su_[fieldi].value(tmVal)/VDash_
);
if (mag(SuValue.value()) <= ROOTVSMALL)
{
// No-op
}
else if (this->useSubMesh())
{
auto tsu = DimensionedField<Type, volMesh>::New
(
name_ + fieldName + "Su",
mesh_,
dimensioned<Type>(SuDims, Zero)
);
UIndirectList<Type>(tsu.ref(), cells_) = SuValue.value();
eqn += tsu;
}
else
{
eqn += SuValue;
}
}
// Implicit source
{
const dimensioned<scalar> SpValue
(
"Sp",
SpDims,
Sp_[fieldi].value(tmVal)/VDash_
);
if (mag(SpValue.value()) <= ROOTVSMALL)
{
// No-op
}
else if (this->useSubMesh())
{
auto tsp = DimensionedField<scalar, volMesh>::New
(
name_ + fieldName + "Sp",
mesh_,
dimensioned<scalar>(SpDims, Zero)
);
UIndirectList<scalar>(tsp.ref(), cells_) = SpValue.value();
eqn += fvm::SuSp(tsp, psi);
}
else
{
eqn += fvm::SuSp(SpValue, psi);
}
}
}
template<class Type>
bool Foam::fv::SemiImplicitSource<Type>::read(const dictionary& dict)
{
VDash_ = 1;
if (fv::cellSetOption::read(dict))
{
volumeMode_ = volumeModeTypeNames_.get("volumeMode", coeffs_);
setFieldData(coeffs_.subDict("injectionRateSuSp"));
// Set volume normalisation
if (volumeMode_ == vmAbsolute)
{
VDash_ = V_;
}
{
setFieldInjectionRates
(
coeffs_.subDict("injectionRateSuSp", keyType::LITERAL)
);
}
return true;
}

86
src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -113,10 +113,9 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef SemiImplicitSource_H
#define SemiImplicitSource_H
#ifndef Foam_SemiImplicitSource_H
#define Foam_SemiImplicitSource_H
#include "Tuple2.H"
#include "cellSetOption.H"
#include "Enum.H"
#include "Function1.H"
@ -128,16 +127,11 @@ namespace Foam
namespace fv
{
// Forward declarations
// Forward Declarations
template<class Type> class SemiImplicitSource;
template<class Type>
Ostream& operator<<
(
Ostream&,
const SemiImplicitSource<Type>&
);
Ostream& operator<<(Ostream&, const SemiImplicitSource<Type>&);
/*---------------------------------------------------------------------------*\
Class SemiImplicitSource Declaration
@ -163,9 +157,9 @@ public:
static const Enum<volumeModeType> volumeModeTypeNames_;
protected:
private:
// Protected Data
// Private Data
//- Volume mode
volumeModeType volumeMode_;
@ -173,15 +167,17 @@ protected:
//- Volume normalisation
scalar VDash_;
//- Source field values
//- Explicit source contributions
PtrList<Function1<Type>> Su_;
//- Linearised implicit contributions
PtrList<Function1<scalar>> Sp_;
// Protected Functions
// Private Member Functions
//- Set the local field data
void setFieldData(const dictionary& dict);
//- Set the source coefficients from "injectionRateSuSp"
void setFieldInjectionRates(const dictionary& dict);
public:
@ -204,40 +200,46 @@ public:
// Member Functions
// Access
// Access
//- Return const access to the volume mode
inline const volumeModeType& volumeMode() const;
//- The volume mode
volumeModeType volumeMode() const noexcept
{
return volumeMode_;
}
// Edit
// Edit
//- Return access to the volume mode
inline volumeModeType& volumeMode();
//- Modifiable access to the volume mode
volumeModeType& volumeMode() noexcept
{
return volumeMode_;
}
// Evaluation
// Evaluation
//- Add explicit contribution to equation
virtual void addSup
(
fvMatrix<Type>& eqn,
const label fieldi
);
//- Add explicit contribution to incompressible equation
virtual void addSup
(
fvMatrix<Type>& eqn,
const label fieldi
);
//- Add explicit contribution to compressible equation
virtual void addSup
(
const volScalarField& rho,
fvMatrix<Type>& eqn,
const label fieldi
);
//- Add explicit contribution to compressible equation
virtual void addSup
(
const volScalarField& rho,
fvMatrix<Type>& eqn,
const label fieldi
);
// IO
// IO
//- Read source dictionary
virtual bool read(const dictionary& dict);
//- Read source dictionary
virtual bool read(const dictionary& dict);
};
@ -254,10 +256,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "SemiImplicitSourceI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

48
src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSourceI.H
View File

@ -1,48 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
-------------------------------------------------------------------------------
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 "SemiImplicitSource.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
inline const typename Foam::fv::SemiImplicitSource<Type>::volumeModeType&
Foam::fv::SemiImplicitSource<Type>::volumeMode() const
{
return volumeMode_;
}
template<class Type>
inline typename Foam::fv::SemiImplicitSource<Type>::volumeModeType&
Foam::fv::SemiImplicitSource<Type>::volumeMode()
{
return volumeMode_;
}
// ************************************************************************* //

0
src/fvOptions/sources/general/semiImplicitSource/semiImplicitSource.C → src/fvOptions/sources/general/semiImplicitSource/semiImplicitSources.C
View File