diff --git a/src/faOptions/faOption/faOptionList.H b/src/faOptions/faOption/faOptionList.H index 18b3b535ab..22a3b04e33 100644 --- a/src/faOptions/faOption/faOptionList.H +++ b/src/faOptions/faOption/faOptionList.H @@ -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; diff --git a/src/finiteArea/faMatrices/faMatricesFwd.H b/src/finiteArea/faMatrices/faMatricesFwd.H index 7383f960d0..15ac5c2bef 100644 --- a/src/finiteArea/faMatrices/faMatricesFwd.H +++ b/src/finiteArea/faMatrices/faMatricesFwd.H @@ -28,8 +28,8 @@ Description \*---------------------------------------------------------------------------*/ -#ifndef faMatricesFwd_H -#define faMatricesFwd_H +#ifndef Foam_faMatricesFwd_H +#define Foam_faMatricesFwd_H #include "fieldTypes.H" diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.C b/src/finiteArea/faMatrices/faMatrix/faMatrix.C index 58e7f8cef6..2307a1bb79 100644 --- a/src/finiteArea/faMatrices/faMatrix/faMatrix.C +++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.C @@ -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::faMatrix << endl; // Initialise coupling coefficients - forAll(psi.mesh().boundary(), patchI) + forAll(psi.mesh().boundary(), patchi) { internalCoeffs_.set ( - patchI, - new Field - ( - psi.mesh().boundary()[patchI].size(), - Zero - ) + patchi, + new Field(psi.mesh().boundary()[patchi].size(), Zero) ); boundaryCoeffs_.set ( - patchI, - new Field - ( - psi.mesh().boundary()[patchI].size(), - Zero - ) + patchi, + new Field(psi.mesh().boundary()[patchi].size(), Zero) ); } @@ -236,7 +228,6 @@ Foam::faMatrix::faMatrix template Foam::faMatrix::faMatrix(const faMatrix& fam) : - refCount(), lduMatrix(fam), psi_(fam.psi_), dimensions_(fam.dimensions_), @@ -246,75 +237,56 @@ Foam::faMatrix::faMatrix(const faMatrix& fam) faceFluxCorrectionPtr_(nullptr) { DebugInFunction - << "copying faMatrix for field " << psi_.name() - << endl; + << "Copying faMatrix for field " << psi_.name() << endl; if (fam.faceFluxCorrectionPtr_) { - faceFluxCorrectionPtr_ = new - GeometricField - ( - *(fam.faceFluxCorrectionPtr_) - ); + faceFluxCorrectionPtr_ = + new GeometricField + ( + *(fam.faceFluxCorrectionPtr_) + ); } } template -Foam::faMatrix::faMatrix -( - const GeometricField& psi, - Istream& is -) +Foam::faMatrix::faMatrix(const tmp>& 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 for field " << psi_.name() - << endl; + << "Copy/Move faMatrix for field " << psi_.name() << endl; - // Initialise coupling coefficients - forAll(psi.mesh().boundary(), patchI) + if (tmat().faceFluxCorrectionPtr_) { - internalCoeffs_.set - ( - patchI, - new Field - ( - psi.mesh().boundary()[patchI].size(), - Zero - ) - ); - - boundaryCoeffs_.set - ( - patchI, - new Field - ( - psi.mesh().boundary()[patchI].size(), - Zero - ) - ); + if (tmat.movable()) + { + faceFluxCorrectionPtr_ = tmat().faceFluxCorrectionPtr_; + tmat().faceFluxCorrectionPtr_ = nullptr; + } + else + { + faceFluxCorrectionPtr_ = + new GeometricField + ( + *(tmat().faceFluxCorrectionPtr_) + ); + } } -} - -template -Foam::tmp> Foam::faMatrix::clone() const -{ - return tmp>::New(*this); + tmat.clear(); } // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * * // - template Foam::faMatrix::~faMatrix() { @@ -601,7 +573,7 @@ void Foam::faMatrix::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::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::operator-=(const tmp>& tfamv) template void Foam::faMatrix::operator+= ( - const GeometricField& su + const DimensionedField& su ) { checkMethod(*this, su, "+="); - source() -= su.mesh().S()*su.internalField(); + source() -= su.mesh().S()*su.field(); +} + + +template +void Foam::faMatrix::operator+= +( + const tmp>& tsu +) +{ + operator+=(tsu()); + tsu.clear(); } @@ -930,11 +913,22 @@ void Foam::faMatrix::operator+= template void Foam::faMatrix::operator-= ( - const GeometricField& su + const DimensionedField& su ) { checkMethod(*this, su, "-="); - source() += su.mesh().S()*su.internalField(); + source() += su.mesh().S()*su.field(); +} + + +template +void Foam::faMatrix::operator-= +( + const tmp>& tsu +) +{ + operator-=(tsu()); + tsu.clear(); } @@ -955,7 +949,7 @@ void Foam::faMatrix::operator+= const dimensioned& su ) { - source() -= su.mesh().S()*su; + source() -= psi().mesh().S()*su; } @@ -965,25 +959,38 @@ void Foam::faMatrix::operator-= const dimensioned& su ) { - source() += su.mesh().S()*su; + source() += psi().mesh().S()*su; } +template +void Foam::faMatrix::operator+=(const Foam::zero) +{} + + +template +void Foam::faMatrix::operator-=(const Foam::zero) +{} + + template void Foam::faMatrix::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::operator*= template void Foam::faMatrix::operator*= ( - const tmp& tvsf + const tmp& tfld ) { - operator*=(tvsf()); - tvsf.clear(); + operator*=(tfld()); + tfld.clear(); +} + + +template +void Foam::faMatrix::operator*= +( + const tmp& tfld +) +{ + operator*=(tfld()); + tfld.clear(); } @@ -1030,34 +1048,32 @@ void Foam::faMatrix::operator*= template void Foam::checkMethod ( - const faMatrix& fam1, - const faMatrix& fam2, + const faMatrix& mat1, + const faMatrix& 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 void Foam::checkMethod ( - const faMatrix& fam, - const GeometricField& vf, + const faMatrix& mat, + const DimensionedField& 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 void Foam::checkMethod ( - const faMatrix& fam, + const faMatrix& mat, const dimensioned& 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 Foam::SolverPerformance Foam::solve ( faMatrix& fam, - Istream& solverControls + const dictionary& solverControls ) { return fam.solve(solverControls); @@ -1127,14 +1141,14 @@ Foam::SolverPerformance Foam::solve template Foam::SolverPerformance Foam::solve ( - const tmp>& tfam, - Istream& solverControls + const tmp>& tmat, + const dictionary& solverControls ) { - SolverPerformance solverPerf = - const_cast&>(tfam()).solve(solverControls); + SolverPerformance solverPerf(tmat.constCast().solve(solverControls)); + + tmat.clear(); - tfam.clear(); return solverPerf; } @@ -1147,12 +1161,12 @@ Foam::SolverPerformance Foam::solve(faMatrix& fam) template -Foam::SolverPerformance Foam::solve(const tmp>& tfam) +Foam::SolverPerformance Foam::solve(const tmp>& tmat) { - SolverPerformance solverPerf = - const_cast&>(tfam()).solve(); + SolverPerformance solverPerf(tmat.constCast().solve()); + + tmat.clear(); - tfam.clear(); return solverPerf; } @@ -1350,12 +1364,12 @@ template Foam::tmp> Foam::operator+ ( const faMatrix& A, - const GeometricField& su + const DimensionedField& su ) { checkMethod(A, su, "+"); - tmp> tC(new faMatrix(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::operator+ template Foam::tmp> Foam::operator+ ( - const tmp>& tA, - const GeometricField& su + const faMatrix& A, + const tmp>& tsu ) { - checkMethod(tA(), su, "+"); - tmp> 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::operator+ ) { checkMethod(A, tsu(), "+"); - tmp> tC(new faMatrix(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 +Foam::tmp> Foam::operator+ +( + const tmp>& tA, + const DimensionedField& su +) +{ + checkMethod(tA(), su, "+"); + tmp> tC(tA.ptr()); + tC.ref().source() -= su.mesh().S()*su.field(); + return tC; +} + + +template +Foam::tmp> Foam::operator+ +( + const tmp>& tA, + const tmp>& tsu +) +{ + checkMethod(tA(), tsu(), "+"); + tmp> tC(tA.ptr()); + tC.ref().source() -= tsu().mesh().S()*tsu().field(); tsu.clear(); return tC; } @@ -1398,7 +1442,7 @@ Foam::tmp> Foam::operator+ { checkMethod(tA(), tsu(), "+"); tmp> 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::operator+ template Foam::tmp> Foam::operator+ ( - const GeometricField& su, + const DimensionedField& su, const faMatrix& A ) { checkMethod(A, su, "+"); - tmp> tC(new faMatrix(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::operator+ template Foam::tmp> Foam::operator+ ( - const GeometricField& su, - const tmp>& tA + const tmp>& tsu, + const faMatrix& A ) { - checkMethod(tA(), su, "+"); - tmp> 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::operator+ ) { checkMethod(A, tsu(), "+"); - tmp> tC(new faMatrix(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 +Foam::tmp> Foam::operator+ +( + const DimensionedField& su, + const tmp>& tA +) +{ + checkMethod(tA(), su, "+"); + tmp> tC(tA.ptr()); + tC.ref().source() -= su.mesh().S()*su.field(); + return tC; +} + + +template +Foam::tmp> Foam::operator+ +( + const tmp>& tsu, + const tmp>& tA +) +{ + checkMethod(tA(), tsu(), "+"); + tmp> tC(tA.ptr()); + tC.ref().source() -= tsu().mesh().S()*tsu().field(); tsu.clear(); return tC; } @@ -1456,7 +1530,7 @@ Foam::tmp> Foam::operator+ { checkMethod(tA(), tsu(), "+"); tmp> 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 Foam::tmp> Foam::operator- ( const faMatrix& A, - const GeometricField& su + const DimensionedField& su ) { checkMethod(A, su, "-"); - tmp> tC(new faMatrix(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::operator- template Foam::tmp> Foam::operator- ( - const tmp>& tA, - const GeometricField& su + const faMatrix& A, + const tmp>& tsu ) { - checkMethod(tA(), su, "-"); - tmp> 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::operator- ) { checkMethod(A, tsu(), "-"); - tmp> tC(new faMatrix(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 +Foam::tmp> Foam::operator- +( + const tmp>& tA, + const DimensionedField& su +) +{ + checkMethod(tA(), su, "-"); + tmp> tC(tA.ptr()); + tC.ref().source() += su.mesh().S()*su.field(); + return tC; +} + + +template +Foam::tmp> Foam::operator- +( + const tmp>& tA, + const tmp>& tsu +) +{ + checkMethod(tA(), tsu(), "-"); + tmp> tC(tA.ptr()); + tC.ref().source() += tsu().mesh().S()*tsu().field(); tsu.clear(); return tC; } @@ -1514,7 +1618,7 @@ Foam::tmp> Foam::operator- { checkMethod(tA(), tsu(), "-"); tmp> 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::operator- template Foam::tmp> Foam::operator- ( - const GeometricField& su, + const DimensionedField& su, const faMatrix& A ) { checkMethod(A, su, "-"); - tmp> tC(new faMatrix(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::operator- template Foam::tmp> Foam::operator- ( - const GeometricField& su, - const tmp>& tA + const tmp>& tsu, + const faMatrix& A ) { - checkMethod(tA(), su, "-"); - tmp> 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::operator- ) { checkMethod(A, tsu(), "-"); - tmp> tC(new faMatrix(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 +Foam::tmp> Foam::operator- +( + const DimensionedField& su, + const tmp>& tA +) +{ + checkMethod(tA(), su, "-"); + tmp> tC(tA.ptr()); + tC.ref().negate(); + tC.ref().source() -= su.mesh().S()*su.field(); + return tC; +} + + +template +Foam::tmp> Foam::operator- +( + const tmp>& tsu, + const tmp>& tA +) +{ + checkMethod(tA(), tsu(), "-"); + tmp> 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::operator- checkMethod(tA(), tsu(), "-"); tmp> 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::operator+ ) { checkMethod(A, su, "+"); - tmp> tC(new faMatrix(A)); + auto tC(A.clone()); tC.ref().source() -= su.value()*A.psi().mesh().S(); return tC; } @@ -1618,7 +1754,7 @@ Foam::tmp> Foam::operator+ ) { checkMethod(A, su, "+"); - tmp> tC(new faMatrix(A)); + auto tC(A.clone()); tC.ref().source() -= su.value()*A.psi().mesh().S(); return tC; } @@ -1646,7 +1782,7 @@ Foam::tmp> Foam::operator- ) { checkMethod(A, su, "-"); - tmp> tC(new faMatrix(A)); + auto tC(A.clone()); tC.ref().source() += su.value()*tC().psi().mesh().S(); return tC; } @@ -1674,7 +1810,7 @@ Foam::tmp> Foam::operator- ) { checkMethod(A, su, "-"); - tmp> tC(new faMatrix(A)); + auto tC(A.clone()); tC.ref().negate(); tC.ref().source() -= su.value()*A.psi().mesh().S(); return tC; @@ -1700,12 +1836,12 @@ template Foam::tmp> Foam::operator== ( const faMatrix& A, - const GeometricField& su + const DimensionedField& su ) { checkMethod(A, su, "=="); - tmp> tC(new faMatrix(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::operator== template Foam::tmp> Foam::operator== ( - const tmp>& tA, - const GeometricField& su + const faMatrix& A, + const tmp>& tsu ) { - checkMethod(tA(), su, "=="); - tmp> 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::operator== ) { checkMethod(A, tsu(), "=="); - tmp> tC(new faMatrix(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 +Foam::tmp> Foam::operator== +( + const tmp>& tA, + const DimensionedField& su +) +{ + checkMethod(tA(), su, "=="); + tmp> tC(tA.ptr()); + tC.ref().source() += su.mesh().S()*su.field(); + return tC; +} + + +template +Foam::tmp> Foam::operator== +( + const tmp>& tA, + const tmp>& tsu +) +{ + checkMethod(tA(), tsu(), "=="); + tmp> tC(tA.ptr()); + tC.ref().source() += tsu().mesh().S()*tsu().field(); tsu.clear(); return tC; } @@ -1748,7 +1914,7 @@ Foam::tmp> Foam::operator== { checkMethod(tA(), tsu(), "=="); tmp> 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::operator== ) { checkMethod(A, su, "=="); - tmp> tC(new faMatrix(A)); + auto tC(A.clone()); tC.ref().source() += A.psi().mesh().S()*su.value(); return tC; } @@ -1782,15 +1948,50 @@ Foam::tmp> Foam::operator== } +template +Foam::tmp> Foam::operator== +( + const faMatrix& A, + const Foam::zero +) +{ + return A; +} + + +template +Foam::tmp> Foam::operator== +( + const tmp>& tA, + const Foam::zero +) +{ + return tA; +} + + template Foam::tmp> Foam::operator* ( - const areaScalarField& vsf, + const areaScalarField::Internal& dsf, const faMatrix& A ) { - tmp> tC(new faMatrix(A)); - tC.ref() *= vsf; + auto tC(A.clone()); + tC.ref() *= dsf; + return tC; +} + + +template +Foam::tmp> Foam::operator* +( + const tmp& tdsf, + const faMatrix& A +) +{ + auto tC(A.clone()); + tC.ref() *= tdsf; return tC; } @@ -1802,7 +2003,7 @@ Foam::tmp> Foam::operator* const faMatrix& A ) { - tmp> tC(new faMatrix(A)); + auto tC(A.clone()); tC.ref() *= tvsf; return tC; } @@ -1811,12 +2012,25 @@ Foam::tmp> Foam::operator* template Foam::tmp> Foam::operator* ( - const areaScalarField& vsf, + const areaScalarField::Internal& dsf, const tmp>& tA ) { tmp> tC(tA.ptr()); - tC.ref() *= vsf; + tC.ref() *= dsf; + return tC; +} + + +template +Foam::tmp> Foam::operator* +( + const tmp& tdsf, + const tmp>& tA +) +{ + tmp> tC(tA.ptr()); + tC.ref() *= tdsf; return tC; } @@ -1841,7 +2055,7 @@ Foam::tmp> Foam::operator* const faMatrix& A ) { - tmp> tC(new faMatrix(A)); + auto tC(A.clone()); tC.ref() *= ds; return tC; } diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.H b/src/finiteArea/faMatrices/faMatrix/faMatrix.H index dcee88c439..b0fc26c2e1 100644 --- a/src/finiteArea/faMatrices/faMatrix/faMatrix.H +++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.H @@ -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 psiFieldType; @@ -88,6 +89,8 @@ public: GeometricField faceFluxFieldType; + /// //- The internal field type for the psi field + /// typedef DimensionedField Internal; private: @@ -219,6 +222,7 @@ public: }; + // Runtime information ClassName("faMatrix"); @@ -234,15 +238,26 @@ public: //- Copy construct faMatrix(const faMatrix&); - //- Construct from Istream given field to solve for + //- Copy/move construct from tmp> + faMatrix(const tmp>&); + + //- 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& psi, Istream& is - ); + ) + : + faMatrix(psi, dimensionSet(is)) + {} - //- Clone - tmp> clone() const; + //- Construct and return a clone + tmp> clone() const + { + return tmp>::New(*this); + } //- Destructor @@ -258,45 +273,45 @@ public: return psi_; } - const dimensionSet& dimensions() const + const dimensionSet& dimensions() const noexcept { return dimensions_; } - Field& source() + Field& source() noexcept { return source_; } - const Field& source() const + const Field& source() const noexcept { return source_; } //- faBoundary scalar field containing pseudo-matrix coeffs //- for internal cells - const FieldField& internalCoeffs() const + const FieldField& internalCoeffs() const noexcept { return internalCoeffs_; } //- faBoundary scalar field containing pseudo-matrix coeffs //- for internal cells - FieldField& internalCoeffs() + FieldField& internalCoeffs() noexcept { return internalCoeffs_; } //- faBoundary scalar field containing pseudo-matrix coeffs //- for boundary cells - const FieldField& boundaryCoeffs() const + const FieldField& boundaryCoeffs() const noexcept { return boundaryCoeffs_; } //- faBoundary scalar field containing pseudo-matrix coeffs //- for boundary cells - FieldField& boundaryCoeffs() + FieldField& 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 solve(const dictionary&); //- Solve returning the solution statistics. @@ -418,6 +433,7 @@ public: void operator=(const faMatrix&); void operator=(const tmp>&); + //- Inplace negate void negate(); void operator+=(const faMatrix&); @@ -426,16 +442,28 @@ public: void operator-=(const faMatrix&); void operator-=(const tmp>&); - void operator+=(const GeometricField&); - void operator+=(const tmp>&); + void operator+=(const DimensionedField&); + void operator+=(const tmp>&); + void operator+= + ( + const tmp>& + ); - void operator-=(const GeometricField&); - void operator-=(const tmp>&); + void operator-=(const DimensionedField&); + void operator-=(const tmp>&); + void operator-= + ( + const tmp>& + ); void operator+=(const dimensioned&); void operator-=(const dimensioned&); - void operator*=(const areaScalarField&); + void operator+=(const Foam::zero); + void operator-=(const Foam::zero); + + void operator*=(const areaScalarField::Internal&); + void operator*=(const tmp&); void operator*=(const tmp&); void operator*=(const dimensioned&); @@ -465,7 +493,7 @@ template void checkMethod ( const faMatrix&, - const GeometricField&, + const DimensionedField&, 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 -SolverPerformance solve(faMatrix&, Istream&); +SolverPerformance solve(faMatrix&, 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 -SolverPerformance solve(const tmp>&, Istream&); +SolverPerformance solve(const tmp>&, const dictionary&); //- Solve returning the solution statistics given convergence tolerance @@ -506,18 +534,21 @@ SolverPerformance solve(const tmp>&); // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // +//- Unary negation template tmp> operator- ( const faMatrix& ); +//- Unary negation template tmp> operator- ( const tmp>& ); + template tmp> operator+ ( @@ -546,6 +577,7 @@ tmp> operator+ const tmp>& ); + template tmp> operator- ( @@ -574,6 +606,7 @@ tmp> operator- const tmp>& ); + template tmp> operator== ( @@ -602,18 +635,19 @@ tmp> operator== const tmp>& ); + template tmp> operator+ ( const faMatrix&, - const GeometricField& + const DimensionedField& ); template tmp> operator+ ( - const tmp>&, - const GeometricField& + const faMatrix&, + const tmp>& ); template @@ -623,6 +657,20 @@ tmp> operator+ const tmp>& ); +template +tmp> operator+ +( + const tmp>&, + const DimensionedField& +); + +template +tmp> operator+ +( + const tmp>&, + const tmp>& +); + template tmp> operator+ ( @@ -633,15 +681,15 @@ tmp> operator+ template tmp> operator+ ( - const GeometricField&, + const DimensionedField&, const faMatrix& ); template tmp> operator+ ( - const GeometricField&, - const tmp>& + const tmp>&, + const faMatrix& ); template @@ -651,6 +699,20 @@ tmp> operator+ const faMatrix& ); +template +tmp> operator+ +( + const DimensionedField&, + const tmp>& +); + +template +tmp> operator+ +( + const tmp>&, + const tmp>& +); + template tmp> operator+ ( @@ -662,14 +724,14 @@ template tmp> operator- ( const faMatrix&, - const GeometricField& + const DimensionedField& ); template tmp> operator- ( - const tmp>&, - const GeometricField& + const faMatrix&, + const tmp>& ); template @@ -683,21 +745,36 @@ template tmp> operator- ( const tmp>&, - const tmp>& + const DimensionedField& ); template tmp> operator- ( - const GeometricField&, + const tmp>&, + const tmp>& +); + +template +tmp> operator- +( + const tmp>&, + const tmp>& +); + + +template +tmp> operator- +( + const DimensionedField&, const faMatrix& ); template tmp> operator- ( - const GeometricField&, - const tmp>& + const tmp>&, + const faMatrix& ); template @@ -707,6 +784,20 @@ tmp> operator- const faMatrix& ); +template +tmp> operator- +( + const DimensionedField&, + const tmp>& +); + +template +tmp> operator- +( + const tmp>&, + const tmp>& +); + template tmp> operator- ( @@ -774,14 +865,14 @@ template tmp> operator== ( const faMatrix&, - const GeometricField& + const DimensionedField& ); template tmp> operator== ( - const tmp>&, - const GeometricField& + const faMatrix&, + const tmp>& ); template @@ -791,6 +882,20 @@ tmp> operator== const tmp>& ); +template +tmp> operator== +( + const tmp>&, + const DimensionedField& +); + +template +tmp> operator== +( + const tmp>&, + const tmp>& +); + template tmp> operator== ( @@ -813,18 +918,33 @@ tmp> operator== ); +template +tmp> operator== +( + const faMatrix&, + const Foam::zero +); + +template +tmp> operator== +( + const tmp>&, + const Foam::zero +); + + template tmp> operator* ( - const areaScalarField&, + const areaScalarField::Internal&, const faMatrix& ); template tmp> operator* ( - const areaScalarField&, - const tmp>& + const tmp&, + const faMatrix& ); template @@ -834,6 +954,20 @@ tmp> operator* const faMatrix& ); +template +tmp> operator* +( + const areaScalarField::Internal&, + const tmp>& +); + +template +tmp> operator* +( + const tmp&, + const tmp>& +); + template tmp> operator* ( @@ -841,6 +975,7 @@ tmp> operator* const tmp>& ); + template tmp> operator* ( diff --git a/src/finiteArea/finiteArea/fam/famSup.C b/src/finiteArea/finiteArea/fam/famSup.C index b0ad234b9a..96747d1b5e 100644 --- a/src/finiteArea/finiteArea/fam/famSup.C +++ b/src/finiteArea/finiteArea/fam/famSup.C @@ -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 -tmp> -Su +Foam::zeroField +Foam::fam::Su ( - const GeometricField& su, - const GeometricField& vf + const Foam::zero, + const GeometricField& fld ) { - const faMesh& mesh = vf.mesh(); - - tmp> tfam - ( - new faMatrix - ( - vf, - dimArea*su.dimensions() - ) - ); - faMatrix& fam = tfam.ref(); - - fam.source() -= mesh.S()*su.internalField(); - - return tfam; + return zeroField(); } + template -tmp> -Su +Foam::tmp> +Foam::fam::Su +( + const dimensioned& su, + const GeometricField& fld +) +{ + auto tmat = tmp>::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 +Foam::tmp> +Foam::fam::Su +( + const DimensionedField& su, + const GeometricField& fld +) +{ + auto tmat = tmp>::New + ( + fld, + dimArea*su.dimensions() + ); + auto& mat = tmat.ref(); + const auto& domain = fld.mesh().S(); + + mat.source() -= domain*su.field(); + + return tmat; +} + + +template +Foam::tmp> +Foam::fam::Su +( + const tmp>& tsu, + const GeometricField& fld +) +{ + tmp> tmat = fam::Su(tsu(), fld); + tsu.clear(); + return tmat; +} + + +template +Foam::tmp> +Foam::fam::Su ( const tmp>& tsu, - const GeometricField& vf + const GeometricField& fld ) { - tmp> tfam = fam::Su(tsu(), vf); + tmp> tmat = fam::Su(tsu(), fld); tsu.clear(); - return tfam; + return tmat; } +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + template -tmp> -Sp +Foam::zeroField +Foam::fam::Sp ( - const areaScalarField& sp, - const GeometricField& vf + const Foam::zero, + const GeometricField& fld ) { - const faMesh& mesh = vf.mesh(); - - tmp> tfam - ( - new faMatrix - ( - vf, - dimArea*sp.dimensions()*vf.dimensions() - ) - ); - faMatrix& fam = tfam.ref(); - - fam.diag() += mesh.S()*sp.internalField(); - - return tfam; -} - -template -tmp> -Sp -( - const tmp& tsp, - const GeometricField& vf -) -{ - tmp> tfam = fam::Sp(tsp(), vf); - tsp.clear(); - return tfam; + return zeroField(); } template -tmp> -Sp +Foam::tmp> +Foam::fam::Sp ( const dimensionedScalar& sp, - const GeometricField& vf + const GeometricField& fld ) { - const faMesh& mesh = vf.mesh(); - - tmp> tfam + auto tmat = tmp>::New ( - new faMatrix - ( - vf, - dimArea*sp.dimensions()*vf.dimensions() - ) + fld, + dimArea*sp.dimensions()*fld.dimensions() ); - faMatrix& 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 -tmp> -SuSp +Foam::tmp> +Foam::fam::Sp ( - const areaScalarField& sp, - const GeometricField& vf + const DimensionedField& sp, + const GeometricField& fld ) { - const faMesh& mesh = vf.mesh(); - - tmp> tfam + auto tmat = tmp>::New ( - new faMatrix - ( - vf, - dimArea*sp.dimensions()*vf.dimensions() - ) + fld, + dimArea*sp.dimensions()*fld.dimensions() ); - faMatrix& 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 -tmp> -SuSp +Foam::tmp> +Foam::fam::Sp +( + const tmp>& tsp, + const GeometricField& fld +) +{ + tmp> tmat = fam::Sp(tsp(), fld); + tsp.clear(); + return tmat; +} + + +template +Foam::tmp> +Foam::fam::Sp ( const tmp& tsp, - const GeometricField& vf + const GeometricField& fld ) { - tmp> tfam = fam::SuSp(tsp(), vf); + tmp> tmat = fam::Sp(tsp(), fld); tsp.clear(); - return tfam; + return tmat; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -} // End namespace fam +template +Foam::zeroField +Foam::fam::SuSp +( + const Foam::zero, + const GeometricField& fld +) +{ + return zeroField(); +} -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -} // End namespace Foam +template +Foam::tmp> +Foam::fam::SuSp +( + const dimensionedScalar& susp, + const GeometricField& fld +) +{ + auto tmat = tmp>::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 +Foam::tmp> +Foam::fam::SuSp +( + const DimensionedField& susp, + const GeometricField& fld +) +{ + auto tmat = tmp>::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 +Foam::tmp> +Foam::fam::SuSp +( + const tmp>& tsusp, + const GeometricField& fld +) +{ + tmp> tmat = fam::SuSp(tsusp(), fld); + tsusp.clear(); + return tmat; +} + + +template +Foam::tmp> +Foam::fam::SuSp +( + const tmp& tsusp, + const GeometricField& fld +) +{ + tmp> tmat = fam::SuSp(tsusp(), fld); + tsusp.clear(); + return tmat; +} + // ************************************************************************* // diff --git a/src/finiteArea/finiteArea/fam/famSup.H b/src/finiteArea/finiteArea/fam/famSup.H index 6779bc0bff..ccc04d98fb 100644 --- a/src/finiteArea/finiteArea/fam/famSup.H +++ b/src/finiteArea/finiteArea/fam/famSup.H @@ -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 + zeroField Su + ( + const Foam::zero, + const GeometricField& + ); + + //- A uniform source (no-op for small values) template tmp> Su ( - const GeometricField&, + const dimensioned& su, + const GeometricField& + ); + + template + tmp> Su + ( + const DimensionedField& su, + const GeometricField& + ); + + template + tmp> Su + ( + const tmp>&, const GeometricField& ); @@ -70,35 +95,73 @@ namespace fam // Implicit source + //- A no-op source + template + zeroField Sp + ( + const Foam::zero, + const GeometricField& + ); + + //- A uniform source (no-op for small values) template tmp> Sp ( - const areaScalarField&, + const dimensionedScalar& sp, const GeometricField& ); template tmp> Sp ( - const tmp&, + const DimensionedField& sp, const GeometricField& ); - template tmp> Sp ( - const dimensionedScalar&, + const tmp>&, + const GeometricField& + ); + + template + tmp> Sp + ( + const tmp>&, const GeometricField& ); // Implicit/Explicit source depending on sign of coefficient + //- A no-op source + template + zeroField SuSp + ( + const Foam::zero, + const GeometricField& + ); + + //- A uniform source (no-op for small values) template tmp> SuSp ( - const areaScalarField&, + const dimensionedScalar& susp, + const GeometricField& + ); + + template + tmp> SuSp + ( + const DimensionedField& susp, + const GeometricField& + ); + + template + tmp> SuSp + ( + const tmp>&, const GeometricField& ); diff --git a/src/finiteVolume/cfdTools/general/fvOptions/fvOptionList.H b/src/finiteVolume/cfdTools/general/fvOptions/fvOptionList.H index 14b8ea08c9..d34f456fbe 100644 --- a/src/finiteVolume/cfdTools/general/fvOptions/fvOptionList.H +++ b/src/finiteVolume/cfdTools/general/fvOptions/fvOptionList.H @@ -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_; diff --git a/src/finiteVolume/finiteVolume/fvm/fvmSup.C b/src/finiteVolume/finiteVolume/fvm/fvmSup.C index a7ce415ff5..b51dbc810b 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmSup.C +++ b/src/finiteVolume/finiteVolume/fvm/fvmSup.C @@ -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 -Foam::tmp> -Foam::fvm::Su -( - const DimensionedField& su, - const GeometricField& vf -) -{ - const fvMesh& mesh = vf.mesh(); - - tmp> tfvm - ( - new fvMatrix - ( - vf, - dimVol*su.dimensions() - ) - ); - fvMatrix& fvm = tfvm.ref(); - - fvm.source() -= mesh.V()*su.field(); - - return tfvm; -} - - -template -Foam::tmp> -Foam::fvm::Su -( - const tmp>& tsu, - const GeometricField& vf -) -{ - tmp> tfvm = fvm::Su(tsu(), vf); - tsu.clear(); - return tfvm; -} - - -template -Foam::tmp> -Foam::fvm::Su -( - const tmp>& tsu, - const GeometricField& vf -) -{ - tmp> tfvm = fvm::Su(tsu(), vf); - tsu.clear(); - return tfvm; -} - - template Foam::zeroField Foam::fvm::Su ( - const zero&, - const GeometricField& vf + const Foam::zero, + const GeometricField& fld ) { return zeroField(); @@ -99,89 +44,86 @@ Foam::fvm::Su template Foam::tmp> -Foam::fvm::Sp +Foam::fvm::Su ( - const volScalarField::Internal& sp, - const GeometricField& vf + const dimensioned& su, + const GeometricField& fld ) { - const fvMesh& mesh = vf.mesh(); - - tmp> tfvm + auto tmat = tmp>::New ( - new fvMatrix - ( - vf, - dimVol*sp.dimensions()*vf.dimensions() - ) + fld, + dimVol*su.dimensions() ); - fvMatrix& 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 Foam::tmp> -Foam::fvm::Sp +Foam::fvm::Su ( - const tmp& tsp, - const GeometricField& vf + const DimensionedField& su, + const GeometricField& fld ) { - tmp> tfvm = fvm::Sp(tsp(), vf); - tsp.clear(); - return tfvm; -} - - -template -Foam::tmp> -Foam::fvm::Sp -( - const tmp& tsp, - const GeometricField& vf -) -{ - tmp> tfvm = fvm::Sp(tsp(), vf); - tsp.clear(); - return tfvm; -} - - -template -Foam::tmp> -Foam::fvm::Sp -( - const dimensionedScalar& sp, - const GeometricField& vf -) -{ - const fvMesh& mesh = vf.mesh(); - - tmp> tfvm + auto tmat = tmp>::New ( - new fvMatrix - ( - vf, - dimVol*sp.dimensions()*vf.dimensions() - ) + fld, + dimVol*su.dimensions() ); - fvMatrix& 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 +Foam::tmp> +Foam::fvm::Su +( + const tmp>& tsu, + const GeometricField& fld +) +{ + tmp> tmat = fvm::Su(tsu(), fld); + tsu.clear(); + return tmat; +} + + +template +Foam::tmp> +Foam::fvm::Su +( + const tmp>& tsu, + const GeometricField& fld +) +{ + tmp> tmat = fvm::Su(tsu(), fld); + tsu.clear(); + return tmat; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + template Foam::zeroField Foam::fvm::Sp ( - const zero&, + const Foam::zero, const GeometricField& ) { @@ -191,30 +133,90 @@ Foam::fvm::Sp template Foam::tmp> -Foam::fvm::SuSp +Foam::fvm::Sp ( - const volScalarField::Internal& susp, - const GeometricField& vf + const dimensionedScalar& sp, + const GeometricField& fld ) { - const fvMesh& mesh = vf.mesh(); - - tmp> tfvm + auto tmat = tmp>::New ( - new fvMatrix - ( - vf, - dimVol*susp.dimensions()*vf.dimensions() - ) + fld, + dimVol*sp.dimensions()*fld.dimensions() ); - fvMatrix& 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 +Foam::tmp> +Foam::fvm::Sp +( + const DimensionedField& sp, + const GeometricField& fld +) +{ + auto tmat = tmp>::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 +Foam::tmp> +Foam::fvm::Sp +( + const tmp>& tsp, + const GeometricField& fld +) +{ + tmp> tmat = fvm::Sp(tsp(), fld); + tsp.clear(); + return tmat; +} + + +template +Foam::tmp> +Foam::fvm::Sp +( + const tmp& tsp, + const GeometricField& fld +) +{ + tmp> tmat = fvm::Sp(tsp(), fld); + tsp.clear(); + return tmat; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +Foam::zeroField +Foam::fvm::SuSp +( + const Foam::zero, + const GeometricField& fld +) +{ + return zeroField(); } @@ -222,13 +224,66 @@ template Foam::tmp> Foam::fvm::SuSp ( - const tmp& tsusp, - const GeometricField& vf + const dimensionedScalar& susp, + const GeometricField& fld ) { - tmp> tfvm = fvm::SuSp(tsusp(), vf); + auto tmat = tmp>::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 +Foam::tmp> +Foam::fvm::SuSp +( + const DimensionedField& susp, + const GeometricField& fld +) +{ + auto tmat = tmp>::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 +Foam::tmp> +Foam::fvm::SuSp +( + const tmp>& tsusp, + const GeometricField& fld +) +{ + tmp> tmat = fvm::SuSp(tsusp(), fld); tsusp.clear(); - return tfvm; + return tmat; } @@ -237,24 +292,12 @@ Foam::tmp> Foam::fvm::SuSp ( const tmp& tsusp, - const GeometricField& vf + const GeometricField& fld ) { - tmp> tfvm = fvm::SuSp(tsusp(), vf); + tmp> tmat = fvm::SuSp(tsusp(), fld); tsusp.clear(); - return tfvm; -} - - -template -Foam::zeroField -Foam::fvm::SuSp -( - const zero&, - const GeometricField& vf -) -{ - return zeroField(); + return tmat; } diff --git a/src/finiteVolume/finiteVolume/fvm/fvmSup.H b/src/finiteVolume/finiteVolume/fvm/fvmSup.H index 8094c8f7c0..45c091fa6c 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmSup.H +++ b/src/finiteVolume/finiteVolume/fvm/fvmSup.H @@ -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 + zeroField Su + ( + const Foam::zero, + const GeometricField& + ); + + //- A uniform source (no-op for small values) template tmp> Su ( - const DimensionedField&, + const dimensioned& su, + const GeometricField& + ); + + template + tmp> Su + ( + const DimensionedField& su, const GeometricField& ); @@ -75,27 +92,36 @@ namespace fvm const GeometricField& ); - template - zeroField Su - ( - const zero&, - const GeometricField& - ); - // Implicit source + //- A no-op source + template + zeroField Sp + ( + const Foam::zero, + const GeometricField& + ); + + //- A uniform source (no-op for small values) template tmp> Sp ( - const volScalarField::Internal&, + const dimensionedScalar& sp, const GeometricField& ); template tmp> Sp ( - const tmp&, + const DimensionedField& sp, + const GeometricField& + ); + + template + tmp> Sp + ( + const tmp>&, const GeometricField& ); @@ -107,35 +133,35 @@ namespace fvm ); - template - tmp> Sp - ( - const dimensionedScalar&, - const GeometricField& - ); - - - template - zeroField Sp - ( - const zero&, - const GeometricField& - ); - - // Implicit/Explicit source depending on sign of coefficient + //- A no-op source + template + zeroField SuSp + ( + const Foam::zero, + const GeometricField& + ); + + //- A uniform source (no-op for small values) template tmp> SuSp ( - const volScalarField::Internal&, + const dimensionedScalar& susp, const GeometricField& ); template tmp> SuSp ( - const tmp&, + const DimensionedField& susp, + const GeometricField& + ); + + template + tmp> SuSp + ( + const tmp>&, const GeometricField& ); @@ -145,13 +171,6 @@ namespace fvm const tmp&, const GeometricField& ); - - template - zeroField SuSp - ( - const zero&, - const GeometricField& - ); } diff --git a/src/finiteVolume/fvMatrices/fvMatricesFwd.H b/src/finiteVolume/fvMatrices/fvMatricesFwd.H index f5a8733c64..7c1a616e32 100644 --- a/src/finiteVolume/fvMatrices/fvMatricesFwd.H +++ b/src/finiteVolume/fvMatrices/fvMatricesFwd.H @@ -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 fvMatrix; +template class fvMatrix; typedef fvMatrix fvScalarMatrix; typedef fvMatrix fvVectorMatrix; diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C index bcefd313e3..1184146fb1 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C @@ -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::checkImplicit(const label fieldi) // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - template Foam::fvMatrix::fvMatrix ( @@ -350,7 +349,7 @@ Foam::fvMatrix::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::fvMatrix internalCoeffs_.set ( patchi, - new Field - ( - psi.mesh().boundary()[patchi].size(), - Zero - ) + new Field(psi.mesh().boundary()[patchi].size(), Zero) ); boundaryCoeffs_.set ( patchi, - new Field - ( - psi.mesh().boundary()[patchi].size(), - Zero - ) + new Field(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::fvMatrix template Foam::fvMatrix::fvMatrix(const fvMatrix& fvm) : - refCount(), lduMatrix(fvm), psi_(fvm.psi_), useImplicit_(fvm.useImplicit_), @@ -423,113 +413,40 @@ Foam::fvMatrix::fvMatrix(const fvMatrix& fvm) template -Foam::fvMatrix::fvMatrix(const tmp>& tfvm) +Foam::fvMatrix::fvMatrix(const tmp>& tmat) : - lduMatrix - ( - const_cast&>(tfvm()), - tfvm.isTmp() - ), - psi_(tfvm().psi_), - useImplicit_(tfvm().useImplicit_), - lduAssemblyName_(tfvm().lduAssemblyName_), - nMatrix_(tfvm().nMatrix_), - dimensions_(tfvm().dimensions_), - source_ - ( - const_cast&>(tfvm()).source_, - tfvm.isTmp() - ), - internalCoeffs_ - ( - const_cast&>(tfvm()).internalCoeffs_, - tfvm.isTmp() - ), - boundaryCoeffs_ - ( - const_cast&>(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 for field " << psi_.name() << endl; + << "Copy/move fvMatrix 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 ( - *(tfvm().faceFluxCorrectionPtr_) + *(tmat().faceFluxCorrectionPtr_) ); } } - tfvm.clear(); -} - - -template -Foam::fvMatrix::fvMatrix -( - const GeometricField& 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 for field " << psi_.name() << endl; - - checkImplicit(); - - // Initialise coupling coefficients - forAll(psi.mesh().boundary(), patchi) - { - internalCoeffs_.set - ( - patchi, - new Field - ( - psi.mesh().boundary()[patchi].size(), - Zero - ) - ); - - boundaryCoeffs_.set - ( - patchi, - new Field - ( - psi.mesh().boundary()[patchi].size(), - Zero - ) - ); - } -} - - -template -Foam::tmp> Foam::fvMatrix::clone() const -{ - return tmp>::New(*this); + tmat.clear(); } @@ -769,21 +686,13 @@ void Foam::fvMatrix::setBounAndInterCoeffs() internalCoeffs_.set ( interfaceI, - new Field - ( - psi.mesh().boundary()[patchi].size(), - Zero - ) + new Field(psi.mesh().boundary()[patchi].size(), Zero) ); boundaryCoeffs_.set ( interfaceI, - new Field - ( - psi.mesh().boundary()[patchi].size(), - Zero - ) + new Field(psi.mesh().boundary()[patchi].size(), Zero) ); interfaceI++; } @@ -1882,18 +1791,12 @@ void Foam::fvMatrix::operator-= template -void Foam::fvMatrix::operator+= -( - const zero& -) +void Foam::fvMatrix::operator+=(const Foam::zero) {} template -void Foam::fvMatrix::operator-= -( - const zero& -) +void Foam::fvMatrix::operator-=(const Foam::zero) {} @@ -1930,22 +1833,22 @@ void Foam::fvMatrix::operator*= template void Foam::fvMatrix::operator*= ( - const tmp& tdsf + const tmp& tfld ) { - operator*=(tdsf()); - tdsf.clear(); + operator*=(tfld()); + tfld.clear(); } template void Foam::fvMatrix::operator*= ( - const tmp& tvsf + const tmp& tfld ) { - operator*=(tvsf()); - tvsf.clear(); + operator*=(tfld()); + tfld.clear(); } @@ -1973,34 +1876,32 @@ void Foam::fvMatrix::operator*= template void Foam::checkMethod ( - const fvMatrix& fvm1, - const fvMatrix& fvm2, + const fvMatrix& mat1, + const fvMatrix& 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 void Foam::checkMethod ( - const fvMatrix& fvm, - const DimensionedField& df, + const fvMatrix& mat, + const DimensionedField& 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 void Foam::checkMethod ( - const fvMatrix& fvm, + const fvMatrix& mat, const dimensioned& 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 Foam::solve template Foam::SolverPerformance Foam::solve ( - const tmp>& tfvm, + const tmp>& tmat, const dictionary& solverControls ) { - SolverPerformance solverPerf = - const_cast&>(tfvm()).solve(solverControls); + SolverPerformance solverPerf(tmat.constCast().solve(solverControls)); - tfvm.clear(); + tmat.clear(); return solverPerf; } @@ -2088,12 +1987,11 @@ Foam::SolverPerformance Foam::solve(fvMatrix& fvm) } template -Foam::SolverPerformance Foam::solve(const tmp>& tfvm) +Foam::SolverPerformance Foam::solve(const tmp>& tmat) { - SolverPerformance solverPerf = - const_cast&>(tfvm()).solve(); + SolverPerformance solverPerf(tmat.constCast().solve()); - tfvm.clear(); + tmat.clear(); return solverPerf; } @@ -2289,7 +2187,7 @@ template Foam::tmp> Foam::operator== ( const fvMatrix& A, - const zero& + const Foam::zero ) { return A; @@ -2300,7 +2198,7 @@ template Foam::tmp> Foam::operator== ( const tmp>& tA, - const zero& + const Foam::zero ) { return tA; diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H index 506ced0edf..c6021a514d 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H @@ -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> operator& ); - /*---------------------------------------------------------------------------*\ Class fvMatrix Declaration \*---------------------------------------------------------------------------*/ @@ -125,7 +124,7 @@ public: // Public Types - //- Field type for psi + //- The geometric field type for psi typedef GeometricField psiFieldType; @@ -135,6 +134,10 @@ public: GeometricField faceFluxFieldType; + /// //- The internal field type for the psi field + /// typedef DimensionedField Internal; + + private: // Private Data @@ -286,6 +289,7 @@ public: }; + // Runtime information ClassName("fvMatrix"); @@ -304,15 +308,23 @@ public: //- Copy/move construct from tmp> fvMatrix(const tmp>&); - //- 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& psi, Istream& is - ); + ) + : + fvMatrix(psi, dimensionSet(is)) + {} - //- Clone - tmp> clone() const; + //- Construct and return a clone + tmp> clone() const + { + return tmp>::New(*this); + } //- Destructor @@ -439,45 +451,45 @@ public: } - const dimensionSet& dimensions() const + const dimensionSet& dimensions() const noexcept { return dimensions_; } - Field& source() + Field& source() noexcept { return source_; } - const Field& source() const + const Field& source() const noexcept { return source_; } //- fvBoundary scalar field containing pseudo-matrix coeffs //- for internal cells - const FieldField& internalCoeffs() const + const FieldField& internalCoeffs() const noexcept { return internalCoeffs_; } //- fvBoundary scalar field containing pseudo-matrix coeffs //- for internal cells - FieldField& internalCoeffs() + FieldField& internalCoeffs() noexcept { return internalCoeffs_; } //- fvBoundary scalar field containing pseudo-matrix coeffs //- for boundary cells - const FieldField& boundaryCoeffs() const + const FieldField& boundaryCoeffs() const noexcept { return boundaryCoeffs_; } //- fvBoundary scalar field containing pseudo-matrix coeffs //- for boundary cells - FieldField& boundaryCoeffs() + FieldField& boundaryCoeffs() noexcept { return boundaryCoeffs_; } @@ -639,6 +651,7 @@ public: void operator=(const fvMatrix&); void operator=(const tmp>&); + //- Inplace negate void negate(); void operator+=(const fvMatrix&); @@ -647,27 +660,15 @@ public: void operator-=(const fvMatrix&); void operator-=(const tmp>&); - void operator+= - ( - const DimensionedField& - ); - void operator+= - ( - const tmp>& - ); + void operator+=(const DimensionedField&); + void operator+=(const tmp>&); void operator+= ( const tmp>& ); - void operator-= - ( - const DimensionedField& - ); - void operator-= - ( - const tmp>& - ); + void operator-=(const DimensionedField&); + void operator-=(const tmp>&); void operator-= ( const tmp>& @@ -676,8 +677,8 @@ public: void operator+=(const dimensioned&); void operator-=(const dimensioned&); - 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&); @@ -888,23 +889,25 @@ template tmp> operator== ( const fvMatrix&, - const zero& + const Foam::zero ); template tmp> operator== ( const tmp>&, - const zero& + const Foam::zero ); +//- Unary negation template tmp> operator- ( const fvMatrix& ); +//- Unary negation template tmp> operator- ( diff --git a/src/fvOptions/Make/files b/src/fvOptions/Make/files index 01e2288baa..534e9cd2d7 100644 --- a/src/fvOptions/Make/files +++ b/src/fvOptions/Make/files @@ -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 diff --git a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C index 777d944635..45698769a7 100644 --- a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C +++ b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C @@ -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::volumeModeTypeNames_ }); -// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template -void Foam::fv::SemiImplicitSource::setFieldData(const dictionary& dict) +void Foam::fv::SemiImplicitSource::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::New("Su", subdict, &mesh_)); + Sp_.set(count, Function1::New("Sp", subdict, &mesh_)); + } + else { Tuple2 injectionRate; dEntry.readEntry(injectionRate); @@ -88,21 +99,10 @@ void Foam::fv::SemiImplicitSource::setFieldData(const dictionary& dict) ) ); } - else - { - const dictionary& Sdict = dEntry.dict(); - Su_.set(count, Function1::New("Su", Sdict, &mesh_)); - Sp_.set(count, Function1::New("Sp", Sdict, &mesh_)); - } + fieldNames_[count] = fieldName; ++count; } - - // Set volume normalisation - if (volumeMode_ == vmAbsolute) - { - VDash_ = V_; - } } @@ -119,7 +119,7 @@ Foam::fv::SemiImplicitSource::SemiImplicitSource : fv::cellSetOption(name, modelType, dict, mesh), volumeMode_(vmAbsolute), - VDash_(1.0) + VDash_(1) { read(dict); } @@ -134,51 +134,7 @@ void Foam::fv::SemiImplicitSource::addSup const label fieldi ) { - if (debug) - { - Info<< "SemiImplicitSource<" << pTraits::typeName - << ">::addSup for source " << name_ << endl; - } - - const GeometricField& psi = eqn.psi(); - - typename GeometricField::Internal Su - ( - IOobject - ( - name_ + fieldNames_[fieldi] + "Su", - mesh_.time().timeName(), - mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh_, - dimensioned(eqn.dimensions()/dimVolume, Zero), - false - ); - - const scalar tmVal = mesh_.time().timeOutputValue(); - - UIndirectList(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(Su.dimensions()/psi.dimensions(), Zero), - false - ); - - UIndirectList(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::addSup << ">::addSup for source " << name_ << endl; } - return this->addSup(eqn, fieldi); + const GeometricField& 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 SuValue + ( + "Su", + SuDims, + Su_[fieldi].value(tmVal)/VDash_ + ); + + if (mag(SuValue.value()) <= ROOTVSMALL) + { + // No-op + } + else if (this->useSubMesh()) + { + auto tsu = DimensionedField::New + ( + name_ + fieldName + "Su", + mesh_, + dimensioned(SuDims, Zero) + ); + UIndirectList(tsu.ref(), cells_) = SuValue.value(); + + eqn += tsu; + } + else + { + eqn += SuValue; + } + } + + + // Implicit source + { + const dimensioned SpValue + ( + "Sp", + SpDims, + Sp_[fieldi].value(tmVal)/VDash_ + ); + + if (mag(SpValue.value()) <= ROOTVSMALL) + { + // No-op + } + else if (this->useSubMesh()) + { + auto tsp = DimensionedField::New + ( + name_ + fieldName + "Sp", + mesh_, + dimensioned(SpDims, Zero) + ); + UIndirectList(tsp.ref(), cells_) = SpValue.value(); + + eqn += fvm::SuSp(tsp, psi); + } + else + { + eqn += fvm::SuSp(SpValue, psi); + } + } } template bool Foam::fv::SemiImplicitSource::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; } diff --git a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H index d3a375dcc4..f51888d00c 100644 --- a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H +++ b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H @@ -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 SemiImplicitSource; template -Ostream& operator<< -( - Ostream&, - const SemiImplicitSource& -); - +Ostream& operator<<(Ostream&, const SemiImplicitSource&); /*---------------------------------------------------------------------------*\ Class SemiImplicitSource Declaration @@ -163,9 +157,9 @@ public: static const Enum 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> Su_; + + //- Linearised implicit contributions PtrList> 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& eqn, - const label fieldi - ); + //- Add explicit contribution to incompressible equation + virtual void addSup + ( + fvMatrix& eqn, + const label fieldi + ); - //- Add explicit contribution to compressible equation - virtual void addSup - ( - const volScalarField& rho, - fvMatrix& eqn, - const label fieldi - ); + //- Add explicit contribution to compressible equation + virtual void addSup + ( + const volScalarField& rho, + fvMatrix& 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 // ************************************************************************* // diff --git a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSourceI.H b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSourceI.H deleted file mode 100644 index 7922e86558..0000000000 --- a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSourceI.H +++ /dev/null @@ -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 . - -\*---------------------------------------------------------------------------*/ - -#include "SemiImplicitSource.H" - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -inline const typename Foam::fv::SemiImplicitSource::volumeModeType& -Foam::fv::SemiImplicitSource::volumeMode() const -{ - return volumeMode_; -} - - -template -inline typename Foam::fv::SemiImplicitSource::volumeModeType& -Foam::fv::SemiImplicitSource::volumeMode() -{ - return volumeMode_; -} - - -// ************************************************************************* // diff --git a/src/fvOptions/sources/general/semiImplicitSource/semiImplicitSource.C b/src/fvOptions/sources/general/semiImplicitSource/semiImplicitSources.C similarity index 100% rename from src/fvOptions/sources/general/semiImplicitSource/semiImplicitSource.C rename to src/fvOptions/sources/general/semiImplicitSource/semiImplicitSources.C