Compare commits

...

11 Commits

Author SHA1 Message Date
28fa11e0fc ENH: finiteArea: modernise and clean up code 2022-12-13 13:58:11 +00:00
ac4b5b6e84 ENH: finiteArea: apply modern IOobject traits 2022-12-13 13:14:30 +00:00
ebb1329980 ENH: finiteArea: use std::unique_ptr to better handle pointers 2022-12-13 13:07:29 +00:00
2e3372d554 ENH: regionFaModels: regroup library constituents 2022-12-13 13:07:29 +00:00
884a7fcc25 ENH: regionFaModel: modernise and clean up code 2022-12-13 13:07:25 +00:00
34b7c07d21 ENH: finiteArea: remove various unused fields/methods 2022-12-13 12:37:07 +00:00
c8e543bad8 ENH: finiteArea: improve robustness in code sections vulnerable to math errors
It has been observed that the finite-area framework is prone to numerical
issues when zero-valued edge lenghts, edge/face normals and face areas exist.

To improve exception handling at identified code sections to gracefully
overcome math errors, the problematic entities are lower-bounded by SMALL.
2022-12-13 12:37:07 +00:00
717a4dcb3f ENH: edgeInterpolation: make skewness-correction optimisation compile-time optional 2022-12-13 12:37:07 +00:00
63242edfee BUG: regionFaModels: report min/max for all processors (fixes #2655) 2022-12-13 12:37:00 +00:00
22a6916676 BUG: processorFaPatch: fix the sign of compound assignment 2022-12-13 12:30:22 +00:00
01e799d70d REVERT: leastSquaresFaVectors: pseudo-inverse shows instabilities for 2-D cases 2022-12-13 12:30:22 +00:00
123 changed files with 2896 additions and 3052 deletions

View File

@ -4,7 +4,7 @@
sqrt
(
2*M_PI*sigma*sqr(aMesh.edgeInterpolation::deltaCoeffs())
*aMesh.edgeInterpolation::deltaCoeffs()
*mag(aMesh.edgeInterpolation::deltaCoeffs())
/rhol
)
).value()*runTime.deltaT().value();

View File

@ -208,19 +208,22 @@ void Foam::interfaceTrackingFvMesh::makeUs() const
}
}
UsPtr_ = new areaVectorField
UsPtr_.reset
(
IOobject
new areaVectorField
(
"Us",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
aMesh(),
dimensioned<vector>(dimVelocity, Zero),
patchFieldTypes
IOobject
(
"Us",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
aMesh(),
dimensioned<vector>(dimVelocity, Zero),
patchFieldTypes
)
);
for (const word& patchName : fixedFreeSurfacePatches_)
@ -266,18 +269,21 @@ void Foam::interfaceTrackingFvMesh::makeFsNetPhi() const
<< abort(FatalError);
}
fsNetPhiPtr_ = new areaScalarField
fsNetPhiPtr_.reset
(
IOobject
new areaScalarField
(
"fsNetPhi",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
aMesh(),
dimensionedScalar(dimVelocity*dimArea, Zero)
IOobject
(
"fsNetPhi",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
aMesh(),
dimensionedScalar(dimVelocity*dimArea, Zero)
)
);
}
@ -305,7 +311,8 @@ void Foam::interfaceTrackingFvMesh::makeControlPoints()
if (controlPointsHeader.typeHeaderOk<vectorIOField>())
{
Info<< "Reading control points" << endl;
controlPointsPtr_ =
controlPointsPtr_.reset
(
new vectorIOField
(
IOobject
@ -316,12 +323,14 @@ void Foam::interfaceTrackingFvMesh::makeControlPoints()
IOobject::MUST_READ,
IOobject::AUTO_WRITE
)
);
)
);
}
else
{
Info<< "Creating new control points" << endl;
controlPointsPtr_ =
controlPointsPtr_.reset
(
new vectorIOField
(
IOobject
@ -333,7 +342,8 @@ void Foam::interfaceTrackingFvMesh::makeControlPoints()
IOobject::AUTO_WRITE
),
aMesh().areaCentres().internalField()
);
)
);
initializeControlPointsPosition();
}
@ -352,10 +362,13 @@ void Foam::interfaceTrackingFvMesh::makeMotionPointsMask()
<< abort(FatalError);
}
motionPointsMaskPtr_ = new labelList
motionPointsMaskPtr_.reset
(
mesh().boundaryMesh()[fsPatchIndex()].nPoints(),
1
new labelList
(
mesh().boundaryMesh()[fsPatchIndex()].nPoints(),
1
)
);
// Mark free surface boundary points
@ -414,19 +427,23 @@ void Foam::interfaceTrackingFvMesh::makeDirections()
<< abort(FatalError);
}
pointsDisplacementDirPtr_ =
pointsDisplacementDirPtr_.reset
(
new vectorField
(
mesh().boundaryMesh()[fsPatchIndex()].nPoints(),
Zero
);
)
);
facesDisplacementDirPtr_ =
facesDisplacementDirPtr_.reset
(
new vectorField
(
mesh().boundaryMesh()[fsPatchIndex()].size(),
Zero
);
)
);
if (!normalMotionDir())
{
@ -457,17 +474,20 @@ void Foam::interfaceTrackingFvMesh::makePhis()
<< abort(FatalError);
}
phisPtr_ = new edgeScalarField
phisPtr_.reset
(
IOobject
new edgeScalarField
(
"phis",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
linearEdgeInterpolate(Us()) & aMesh().Le()
IOobject
(
"phis",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
linearEdgeInterpolate(Us()) & aMesh().Le()
)
);
}
@ -484,21 +504,23 @@ void Foam::interfaceTrackingFvMesh::makeSurfactConc() const
<< abort(FatalError);
}
surfactConcPtr_ = new areaScalarField
surfactConcPtr_.reset
(
IOobject
new areaScalarField
(
"Cs",
mesh().time().timeName
IOobject
(
mesh().time().startTime().value()
"Cs",
mesh().time().timeName
(
mesh().time().startTime().value()
),
mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
// mesh().time().timeName(),
mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
aMesh()
aMesh()
)
);
}
@ -515,21 +537,23 @@ void Foam::interfaceTrackingFvMesh::makeBulkSurfactConc() const
<< abort(FatalError);
}
bulkSurfactConcPtr_ = new volScalarField
bulkSurfactConcPtr_.reset
(
IOobject
new volScalarField
(
"C",
mesh().time().timeName
IOobject
(
mesh().time().startTime().value()
"C",
mesh().time().timeName
(
mesh().time().startTime().value()
),
mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
// mesh().time().timeName(),
mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh()
mesh()
)
);
volScalarField& bulkSurfactConc = *bulkSurfactConcPtr_;
@ -554,17 +578,20 @@ void Foam::interfaceTrackingFvMesh::makeSurfaceTension() const
<< abort(FatalError);
}
surfaceTensionPtr_ = new areaScalarField
surfaceTensionPtr_.reset
(
IOobject
new areaScalarField
(
"surfaceTension",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sigma() + surfactant().dSigma(surfactantConcentration())/rho_
IOobject
(
"surfaceTension",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sigma() + surfactant().dSigma(surfactantConcentration())/rho_
)
);
}
@ -584,7 +611,7 @@ void Foam::interfaceTrackingFvMesh::makeSurfactant() const
const dictionary& surfactProp =
motion().subDict("surfactantProperties");
surfactantPtr_ = new surfactantProperties(surfactProp);
surfactantPtr_.reset(new surfactantProperties(surfactProp));
}
@ -613,7 +640,8 @@ void Foam::interfaceTrackingFvMesh::makeContactAngle()
{
Info<< "Reading contact angle field" << endl;
contactAnglePtr_ =
contactAnglePtr_.reset
(
new areaScalarField
(
IOobject
@ -625,7 +653,8 @@ void Foam::interfaceTrackingFvMesh::makeContactAngle()
IOobject::AUTO_WRITE
),
aMesh()
);
)
);
}
}
@ -1643,18 +1672,18 @@ Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh
Foam::interfaceTrackingFvMesh::~interfaceTrackingFvMesh()
{
deleteDemandDrivenData(UsPtr_);
deleteDemandDrivenData(controlPointsPtr_);
deleteDemandDrivenData(motionPointsMaskPtr_);
deleteDemandDrivenData(pointsDisplacementDirPtr_);
deleteDemandDrivenData(facesDisplacementDirPtr_);
deleteDemandDrivenData(fsNetPhiPtr_);
deleteDemandDrivenData(phisPtr_);
deleteDemandDrivenData(surfactConcPtr_);
deleteDemandDrivenData(bulkSurfactConcPtr_);
deleteDemandDrivenData(surfaceTensionPtr_);
deleteDemandDrivenData(surfactantPtr_);
deleteDemandDrivenData(contactAnglePtr_);
UsPtr_.reset(nullptr);
controlPointsPtr_.reset(nullptr);
motionPointsMaskPtr_.reset(nullptr);
pointsDisplacementDirPtr_.reset(nullptr);
facesDisplacementDirPtr_.reset(nullptr);
fsNetPhiPtr_.reset(nullptr);
phisPtr_.reset(nullptr);
surfactConcPtr_.reset(nullptr);
bulkSurfactConcPtr_.reset(nullptr);
surfaceTensionPtr_.reset(nullptr);
surfactantPtr_.reset(nullptr);
contactAnglePtr_.reset(nullptr);
}

View File

@ -107,42 +107,43 @@ class interfaceTrackingFvMesh
label timeIndex_;
//- Free-surface velocity field
mutable areaVectorField* UsPtr_;
mutable std::unique_ptr<areaVectorField> UsPtr_;
//- Points which are attached to the free-surface A side faces
// and which defines the free-surface shape
mutable vectorIOField* controlPointsPtr_;
mutable std::unique_ptr<vectorIOField> controlPointsPtr_;
//- Field which additionally determines
// the motion of free-surface points
mutable labelList* motionPointsMaskPtr_;
mutable std::unique_ptr<labelList> motionPointsMaskPtr_;
//- Displacement direction of free-surface points
mutable vectorField* pointsDisplacementDirPtr_;
mutable std::unique_ptr<vectorField> pointsDisplacementDirPtr_;
//- Displacement direction of free-surface control points
mutable vectorField* facesDisplacementDirPtr_;
mutable std::unique_ptr<vectorField> facesDisplacementDirPtr_;
//- Free-surface net flux
mutable areaScalarField* fsNetPhiPtr_;
mutable std::unique_ptr<areaScalarField> fsNetPhiPtr_;
//- Free-surface flux
mutable edgeScalarField* phisPtr_;
mutable std::unique_ptr<edgeScalarField> phisPtr_;
//- Free-surface surfactant concetration
mutable areaScalarField* surfactConcPtr_;
mutable std::unique_ptr<areaScalarField> surfactConcPtr_;
//- Volume surfactant concetration
mutable volScalarField* bulkSurfactConcPtr_;
mutable std::unique_ptr<volScalarField> bulkSurfactConcPtr_;
//- Surface tension field
mutable areaScalarField* surfaceTensionPtr_;
mutable std::unique_ptr<areaScalarField> surfaceTensionPtr_;
//- Surfactant properties
mutable surfactantProperties* surfactantPtr_;
mutable std::unique_ptr<surfactantProperties> surfactantPtr_;
//- Contact angle
mutable areaScalarField* contactAnglePtr_;
mutable std::unique_ptr<areaScalarField> contactAnglePtr_;
// Private Member Functions
@ -425,7 +426,7 @@ public:
//- Clear control points
void clearControlPoints()
{
deleteDemandDrivenData(controlPointsPtr_);
controlPointsPtr_.reset(nullptr);
}
//- Write VTK freeSurface mesh

View File

@ -49,7 +49,7 @@ void Foam::faMatrix<Type>::addToInternalField
{
FatalErrorInFunction
<< "addressing (" << addr.size()
<< ") and field (" << pf.size() << ") are different sizes" << endl
<< ") and field (" << pf.size() << ") are different sizes"
<< abort(FatalError);
}
@ -87,7 +87,7 @@ void Foam::faMatrix<Type>::subtractFromInternalField
{
FatalErrorInFunction
<< "addressing (" << addr.size()
<< ") and field (" << pf.size() << ") are different sizes" << endl
<< ") and field (" << pf.size() << ") are different sizes"
<< abort(FatalError);
}

View File

@ -234,17 +234,17 @@ void Foam::faMesh::clearGeomNotAreas() const
polyPatchFacesPtr_.reset(nullptr);
polyPatchIdsPtr_.reset(nullptr);
bndConnectPtr_.reset(nullptr);
deleteDemandDrivenData(SPtr_);
deleteDemandDrivenData(patchStartsPtr_);
deleteDemandDrivenData(LePtr_);
deleteDemandDrivenData(magLePtr_);
deleteDemandDrivenData(faceCentresPtr_);
deleteDemandDrivenData(edgeCentresPtr_);
deleteDemandDrivenData(faceAreaNormalsPtr_);
deleteDemandDrivenData(edgeAreaNormalsPtr_);
Sptr_.reset(nullptr);
patchStartsPtr_.reset(nullptr);
LePtr_.reset(nullptr);
magLePtr_.reset(nullptr);
faceCentresPtr_.reset(nullptr);
edgeCentresPtr_.reset(nullptr);
faceAreaNormalsPtr_.reset(nullptr);
edgeAreaNormalsPtr_.reset(nullptr);
pointAreaNormalsPtr_.reset(nullptr);
deleteDemandDrivenData(faceCurvaturesPtr_);
deleteDemandDrivenData(edgeTransformTensorsPtr_);
faceCurvaturesPtr_.reset(nullptr);
edgeTransformTensorsPtr_.reset(nullptr);
}
@ -253,9 +253,9 @@ void Foam::faMesh::clearGeom() const
DebugInFunction << "Clearing geometry" << endl;
clearGeomNotAreas();
deleteDemandDrivenData(S0Ptr_);
deleteDemandDrivenData(S00Ptr_);
deleteDemandDrivenData(correctPatchPointNormalsPtr_);
S0ptr_.reset(nullptr);
S00ptr_.reset(nullptr);
correctPatchPointNormalsPtr_.reset(nullptr);
}
@ -263,7 +263,7 @@ void Foam::faMesh::clearAddressing() const
{
DebugInFunction << "Clearing addressing" << endl;
deleteDemandDrivenData(lduPtr_);
lduPtr_.reset(nullptr);
}
@ -377,9 +377,9 @@ Foam::faMesh::faMesh
bndConnectPtr_(nullptr),
lduPtr_(nullptr),
SPtr_(nullptr),
S0Ptr_(nullptr),
S00Ptr_(nullptr),
Sptr_(nullptr),
S0ptr_(nullptr),
S00ptr_(nullptr),
patchStartsPtr_(nullptr),
LePtr_(nullptr),
magLePtr_(nullptr),
@ -408,19 +408,22 @@ Foam::faMesh::faMesh
if (doInit && fileHandler().isFile(pMesh.time().timePath()/"S0"))
{
S0Ptr_ = new DimensionedField<scalar, areaMesh>
S0ptr_.reset
(
IOobject
new DimensionedField<scalar, areaMesh>
(
"S0",
time().timeName(),
faMesh::meshSubDir,
mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
*this
);
IOobject
(
"S0",
time().timeName(),
faMesh::meshSubDir,
mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
*this
)
);
}
}
@ -483,9 +486,9 @@ Foam::faMesh::faMesh
bndConnectPtr_(nullptr),
lduPtr_(nullptr),
SPtr_(nullptr),
S0Ptr_(nullptr),
S00Ptr_(nullptr),
Sptr_(nullptr),
S0ptr_(nullptr),
S00ptr_(nullptr),
patchStartsPtr_(nullptr),
LePtr_(nullptr),
magLePtr_(nullptr),
@ -564,9 +567,9 @@ Foam::faMesh::faMesh
bndConnectPtr_(nullptr),
lduPtr_(nullptr),
SPtr_(nullptr),
S0Ptr_(nullptr),
S00Ptr_(nullptr),
Sptr_(nullptr),
S0ptr_(nullptr),
S00ptr_(nullptr),
patchStartsPtr_(nullptr),
LePtr_(nullptr),
magLePtr_(nullptr),
@ -655,18 +658,21 @@ Foam::faMesh::faMesh
if (doInit && fileHandler().isFile(pMesh.time().timePath()/"S0"))
{
S0Ptr_ = new DimensionedField<scalar, areaMesh>
S0ptr_.reset
(
IOobject
new DimensionedField<scalar, areaMesh>
(
"S0",
time().timeName(),
faMesh::meshSubDir,
mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
*this
IOobject
(
"S0",
time().timeName(),
faMesh::meshSubDir,
mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
*this
)
);
}
}
@ -813,51 +819,54 @@ const Foam::edgeVectorField& Foam::faMesh::edgeCentres() const
const Foam::DimensionedField<Foam::scalar, Foam::areaMesh>&
Foam::faMesh::S() const
{
if (!SPtr_)
if (!Sptr_)
{
calcS();
}
return *SPtr_;
return *Sptr_;
}
const Foam::DimensionedField<Foam::scalar, Foam::areaMesh>&
Foam::faMesh::S0() const
{
if (!S0Ptr_)
if (!S0ptr_)
{
FatalErrorInFunction
<< "S0 is not available"
<< abort(FatalError);
}
return *S0Ptr_;
return *S0ptr_;
}
const Foam::DimensionedField<Foam::scalar, Foam::areaMesh>&
Foam::faMesh::S00() const
{
if (!S00Ptr_)
if (!S00ptr_)
{
S00Ptr_ = new DimensionedField<scalar, areaMesh>
S00ptr_.reset
(
IOobject
new DimensionedField<scalar, areaMesh>
(
"S00",
time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
S0()
IOobject
(
"S00",
time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
S0()
)
);
S0Ptr_->writeOpt(IOobject::AUTO_WRITE);
S0ptr_->writeOpt(IOobject::AUTO_WRITE);
}
return *S00Ptr_;
return *S00ptr_;
}
@ -954,33 +963,36 @@ bool Foam::faMesh::movePoints()
// Grab old time areas if the time has been incremented
if (curTimeIndex_ < time().timeIndex())
{
if (S00Ptr_ && S0Ptr_)
if (S00ptr_ && S0ptr_)
{
DebugInfo<< "Copy old-old S" << endl;
*S00Ptr_ = *S0Ptr_;
*S00ptr_ = *S0ptr_;
}
if (S0Ptr_)
if (S0ptr_)
{
DebugInfo<< "Copy old S" << endl;
*S0Ptr_ = S();
*S0ptr_ = S();
}
else
{
DebugInfo<< "Creating old cell volumes." << endl;
S0Ptr_ = new DimensionedField<scalar, areaMesh>
S0ptr_.reset
(
IOobject
new DimensionedField<scalar, areaMesh>
(
"S0",
time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
S()
IOobject
(
"S0",
time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
S()
)
);
}
@ -1026,7 +1038,10 @@ Foam::boolList& Foam::faMesh::correctPatchPointNormals() const
{
if (!correctPatchPointNormalsPtr_)
{
correctPatchPointNormalsPtr_ = new boolList(boundary().size(), false);
correctPatchPointNormalsPtr_.reset
(
new boolList(boundary().size(), false)
);
}
return *correctPatchPointNormalsPtr_;

View File

@ -263,52 +263,53 @@ class faMesh
mutable std::unique_ptr<List<labelPair>> bndConnectPtr_;
//- Ldu addressing data
mutable faMeshLduAddressing* lduPtr_;
mutable std::unique_ptr<faMeshLduAddressing> lduPtr_;
// Geometric Data
//- Face areas
mutable DimensionedField<scalar, areaMesh>* SPtr_;
mutable std::unique_ptr<DimensionedField<scalar, areaMesh>> Sptr_;
//- Face areas old time level
mutable DimensionedField<scalar, areaMesh>* S0Ptr_;
mutable std::unique_ptr<DimensionedField<scalar, areaMesh>> S0ptr_;
//- Face areas old-old time level
mutable DimensionedField<scalar, areaMesh>* S00Ptr_;
mutable std::unique_ptr<DimensionedField<scalar, areaMesh>> S00ptr_;
//- Patch starts in the edge list
mutable labelList* patchStartsPtr_;
mutable std::unique_ptr<labelList> patchStartsPtr_;
//- Edge length vectors
mutable edgeVectorField* LePtr_;
mutable std::unique_ptr<edgeVectorField> LePtr_;
//- Mag edge length vectors
mutable edgeScalarField* magLePtr_;
mutable std::unique_ptr<edgeScalarField> magLePtr_;
//- Face centres
mutable areaVectorField* faceCentresPtr_;
mutable std::unique_ptr<areaVectorField> faceCentresPtr_;
//- Edge centres
mutable edgeVectorField* edgeCentresPtr_;
mutable std::unique_ptr<edgeVectorField> edgeCentresPtr_;
//- Face area normals
mutable areaVectorField* faceAreaNormalsPtr_;
mutable std::unique_ptr<areaVectorField> faceAreaNormalsPtr_;
//- Edge area normals
mutable edgeVectorField* edgeAreaNormalsPtr_;
mutable std::unique_ptr<edgeVectorField> edgeAreaNormalsPtr_;
//- Point area normals
mutable std::unique_ptr<vectorField> pointAreaNormalsPtr_;
//- Face curvatures
mutable areaScalarField* faceCurvaturesPtr_;
mutable std::unique_ptr<areaScalarField> faceCurvaturesPtr_;
//- Edge transformation tensors
mutable FieldField<Field, tensor>* edgeTransformTensorsPtr_;
mutable std::unique_ptr<FieldField<Field, tensor>>
edgeTransformTensorsPtr_;
//- Whether point normals must be corrected for a patch
mutable boolList* correctPatchPointNormalsPtr_;
mutable std::unique_ptr<boolList> correctPatchPointNormalsPtr_;
// Other mesh-related data
@ -804,7 +805,7 @@ public:
//- Has face areas: S()
bool hasFaceAreas() const noexcept { return bool(SPtr_); }
bool hasFaceAreas() const noexcept { return bool(Sptr_); }
//- Has face centres: areaCentres()
bool hasAreaCentres() const noexcept { return bool(faceCentresPtr_); }

View File

@ -221,7 +221,7 @@ void Foam::faMesh::calcLduAddressing() const
<< abort(FatalError);
}
lduPtr_ = new faMeshLduAddressing(*this);
lduPtr_.reset(new faMeshLduAddressing(*this));
}
@ -237,7 +237,7 @@ void Foam::faMesh::calcPatchStarts() const
<< abort(FatalError);
}
patchStartsPtr_ = new labelList(boundary().patchStarts());
patchStartsPtr_.reset(new labelList(boundary().patchStarts()));
}
@ -594,6 +594,12 @@ Foam::tmp<Foam::vectorField> Foam::faMesh::calcRawEdgeNormals(int order) const
edgeNormals[edgei].removeCollinear(edgeLine.unitVec());
edgeNormals[edgei].normalise();
// Do not allow any mag(val) < SMALL
if (mag(edgeNormals[edgei]) < SMALL)
{
edgeNormals[edgei] = vector::uniform(SMALL);
}
}
return tedgeNormals;
@ -612,7 +618,8 @@ void Foam::faMesh::calcLe() const
<< abort(FatalError);
}
LePtr_ =
LePtr_.reset
(
new edgeVectorField
(
IOobject
@ -625,8 +632,8 @@ void Foam::faMesh::calcLe() const
*this,
dimLength
// -> calculatedType()
);
)
);
edgeVectorField& Le = *LePtr_;
// Need face centres
@ -658,6 +665,12 @@ void Foam::faMesh::calcLe() const
edges_[edgei].line(localPoints),
edgeNormals[edgei]
);
// Do not allow any mag(val) < SMALL
if (mag(leVectors[edgei]) < SMALL)
{
leVectors[edgei] = vector::uniform(SMALL);
}
}
// Copy internal field
@ -672,6 +685,15 @@ void Foam::faMesh::calcLe() const
{
const faPatch& fap = boundary()[patchi];
bfld[patchi] = fap.patchRawSlice(leVectors);
for (auto& patchEdge : bfld[patchi])
{
// Do not allow any mag(val) < SMALL
if (mag(patchEdge) < SMALL)
{
patchEdge = vector::uniform(SMALL);
}
}
}
}
else
@ -692,6 +714,12 @@ void Foam::faMesh::calcLe() const
edges_[edgei].line(localPoints),
edgeNormals[edgei]
);
// Do not allow any mag(val) < SMALL
if (mag(fld[edgei]) < SMALL)
{
fld[edgei] = vector::uniform(SMALL);
}
}
}
@ -715,6 +743,12 @@ void Foam::faMesh::calcLe() const
bndEdgeNormals[patchEdgei]
);
// Do not allow any mag(val) < SMALL
if (mag(pfld[patchEdgei]) < SMALL)
{
pfld[patchEdgei] = vector::uniform(SMALL);
}
++edgei;
}
}
@ -734,7 +768,8 @@ void Foam::faMesh::calcMagLe() const
<< abort(FatalError);
}
magLePtr_ =
magLePtr_.reset
(
new edgeScalarField
(
IOobject
@ -746,8 +781,8 @@ void Foam::faMesh::calcMagLe() const
),
*this,
dimLength
);
)
);
edgeScalarField& magLe = *magLePtr_;
const pointField& localPoints = points();
@ -759,6 +794,13 @@ void Foam::faMesh::calcMagLe() const
for (const edge& e : internalEdges())
{
*iter = e.mag(localPoints);
// Do not allow any mag(val) < SMALL
if (mag(*iter) < SMALL)
{
*iter = SMALL;
}
++iter;
}
}
@ -774,6 +816,13 @@ void Foam::faMesh::calcMagLe() const
for (const edge& e : boundary()[patchi].patchSlice(edges_))
{
*iter = e.mag(localPoints);
// Do not allow any mag(val) < SMALL
if (mag(*iter) < SMALL)
{
*iter = SMALL;
}
++iter;
}
}
@ -793,7 +842,8 @@ void Foam::faMesh::calcFaceCentres() const
<< abort(FatalError);
}
faceCentresPtr_ =
faceCentresPtr_.reset
(
new areaVectorField
(
IOobject
@ -806,8 +856,8 @@ void Foam::faMesh::calcFaceCentres() const
*this,
dimLength
// -> calculatedType()
);
)
);
areaVectorField& centres = *faceCentresPtr_;
// Need local points
@ -866,7 +916,8 @@ void Foam::faMesh::calcEdgeCentres() const
<< abort(FatalError);
}
edgeCentresPtr_ =
edgeCentresPtr_.reset
(
new edgeVectorField
(
IOobject
@ -879,8 +930,8 @@ void Foam::faMesh::calcEdgeCentres() const
*this,
dimLength
// -> calculatedType()
);
)
);
edgeVectorField& centres = *edgeCentresPtr_;
// Need local points
@ -921,27 +972,30 @@ void Foam::faMesh::calcS() const
DebugInFunction
<< "Calculating areas" << endl;
if (SPtr_)
if (Sptr_)
{
FatalErrorInFunction
<< "S() already allocated"
<< abort(FatalError);
}
SPtr_ = new DimensionedField<scalar, areaMesh>
Sptr_.reset
(
IOobject
new DimensionedField<scalar, areaMesh>
(
"S",
time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
*this,
dimArea
IOobject
(
"S",
time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
*this,
dimArea
)
);
auto& areas = *SPtr_;
auto& areas = *Sptr_;
// No access to fvMesh::magSf(), only polyMesh::faceAreas()
@ -955,6 +1009,12 @@ void Foam::faMesh::calcS() const
forAll(fld, facei)
{
fld[facei] = Foam::mag(meshFaceAreas[facei]);
// Do not allow any mag(val) < SMALL
if (mag(fld[facei]) < SMALL)
{
fld[facei] = SMALL;
}
}
}
else
@ -968,6 +1028,13 @@ void Foam::faMesh::calcS() const
for (const face& f : faces())
{
*iter = f.mag(localPoints);
// Do not allow any mag(val) < SMALL
if (mag(*iter) < SMALL)
{
*iter = SMALL;
}
++iter;
}
}
@ -986,7 +1053,8 @@ void Foam::faMesh::calcFaceAreaNormals() const
<< abort(FatalError);
}
faceAreaNormalsPtr_ =
faceAreaNormalsPtr_.reset
(
new areaVectorField
(
IOobject
@ -998,8 +1066,8 @@ void Foam::faMesh::calcFaceAreaNormals() const
),
*this,
dimless
);
)
);
areaVectorField& faceNormals = *faceAreaNormalsPtr_;
const pointField& localPoints = points();
@ -1028,6 +1096,15 @@ void Foam::faMesh::calcFaceAreaNormals() const
// Make unit normals
fld.normalise();
for (auto& f : fld)
{
// Do not allow any mag(val) < SMALL
if (mag(f) < SMALL)
{
f = vector::uniform(SMALL);
}
}
}
@ -1056,7 +1133,8 @@ void Foam::faMesh::calcEdgeAreaNormals() const
<< abort(FatalError);
}
edgeAreaNormalsPtr_ =
edgeAreaNormalsPtr_.reset
(
new edgeVectorField
(
IOobject
@ -1068,7 +1146,8 @@ void Foam::faMesh::calcEdgeAreaNormals() const
),
*this,
dimless
);
)
);
edgeVectorField& edgeAreaNormals = *edgeAreaNormalsPtr_;
@ -1122,6 +1201,12 @@ void Foam::faMesh::calcEdgeAreaNormals() const
fld[edgei].removeCollinear(edgeLine.unitVec());
fld[edgei].normalise();
// Do not allow any mag(val) < SMALL
if (mag(fld[edgei]) < SMALL)
{
fld[edgei] = vector::uniform(SMALL);
}
}
}
@ -1148,6 +1233,12 @@ void Foam::faMesh::calcEdgeAreaNormals() const
pfld[patchEdgei].removeCollinear(edgeLine.unitVec());
pfld[patchEdgei].normalise();
// Do not allow any mag(val) < SMALL
if (mag(pfld[patchEdgei]) < SMALL)
{
pfld[patchEdgei] = vector::uniform(SMALL);
}
}
}
}
@ -1166,7 +1257,8 @@ void Foam::faMesh::calcFaceCurvatures() const
<< abort(FatalError);
}
faceCurvaturesPtr_ =
faceCurvaturesPtr_.reset
(
new areaScalarField
(
IOobject
@ -1178,8 +1270,8 @@ void Foam::faMesh::calcFaceCurvatures() const
),
*this,
dimless/dimLength
);
)
);
areaScalarField& faceCurvatures = *faceCurvaturesPtr_;
@ -1205,7 +1297,7 @@ void Foam::faMesh::calcEdgeTransformTensors() const
<< abort(FatalError);
}
edgeTransformTensorsPtr_ = new FieldField<Field, tensor>(nEdges_);
edgeTransformTensorsPtr_.reset(new FieldField<Field, tensor>(nEdges_));
auto& edgeTransformTensors = *edgeTransformTensorsPtr_;
// Initialize all transformation tensors
@ -1760,14 +1852,14 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
}
forAll(boundary(), patchI)
forAll(boundary(), patchi)
{
const faPatch& fap = boundary()[patchI];
const faPatch& fap = boundary()[patchi];
if (Pstream::parRun() && isA<processorFaPatch>(fap))
{
const processorFaPatch& procPatch =
refCast<const processorFaPatch>(boundary()[patchI]);
refCast<const processorFaPatch>(boundary()[patchi]);
const labelList& patchPointLabels = procPatch.pointLabels();

View File

@ -61,10 +61,10 @@ void Foam::faAreaMapper::calcAddressing() const
// Prepare a list of new face labels and (preliminary) addressing
// Note: dimensioned to number of boundary faces of polyMesh
newFaceLabelsPtr_ = new labelList(mesh_().nBoundaryFaces(), -1);
newFaceLabelsPtr_.reset(new labelList(mesh_().nBoundaryFaces(), -1));
labelList& newFaceLabels = *newFaceLabelsPtr_;
newFaceLabelsMapPtr_ = new labelList(mesh_().nBoundaryFaces(), -1);
newFaceLabelsMapPtr_.reset(new labelList(mesh_().nBoundaryFaces(), -1));
labelList& newFaceLabelsMap = *newFaceLabelsMapPtr_;
label nNewFaces = 0;
@ -94,7 +94,7 @@ void Foam::faAreaMapper::calcAddressing() const
// Direct mapping: no further faces to add. Resize list
newFaceLabels.setSize(nNewFaces);
directAddrPtr_ = new labelList(newFaceLabels.size());
directAddrPtr_.reset(new labelList(newFaceLabels.size()));
labelList& addr = *directAddrPtr_;
// Adjust for creation of a boundary face from an internal face
@ -114,10 +114,10 @@ void Foam::faAreaMapper::calcAddressing() const
{
// There are further faces to add. Prepare interpolation addressing
// and weights to full size
interpolationAddrPtr_ = new labelListList(newFaceLabels.size());
interpolationAddrPtr_.reset(new labelListList(newFaceLabels.size()));
labelListList& addr = *interpolationAddrPtr_;
weightsPtr_ = new scalarListList(newFaceLabels.size());
weightsPtr_.reset(new scalarListList(newFaceLabels.size()));
scalarListList& w = *weightsPtr_;
// Insert single addressing and weights
@ -256,20 +256,20 @@ void Foam::faAreaMapper::calcAddressing() const
// Inserted objects cannot appear in the new faMesh as they have no master
// HJ, 10/Aug/2011
insertedObjectLabelsPtr_ = new labelList(0);
insertedObjectLabelsPtr_.reset(new labelList(0));
}
void Foam::faAreaMapper::clearOut()
{
deleteDemandDrivenData(newFaceLabelsPtr_);
deleteDemandDrivenData(newFaceLabelsMapPtr_);
newFaceLabelsPtr_.reset(nullptr);
newFaceLabelsMapPtr_.reset(nullptr);
deleteDemandDrivenData(directAddrPtr_);
deleteDemandDrivenData(interpolationAddrPtr_);
deleteDemandDrivenData(weightsPtr_);
directAddrPtr_.reset(nullptr);
interpolationAddrPtr_.reset(nullptr);
weightsPtr_.reset(nullptr);
deleteDemandDrivenData(insertedObjectLabelsPtr_);
insertedObjectLabelsPtr_.reset(nullptr);
hasUnmapped_ = false;
}

View File

@ -80,23 +80,23 @@ class faAreaMapper
label sizeBeforeMapping_;
//- New face labels after mapping
mutable labelList* newFaceLabelsPtr_;
mutable std::unique_ptr<labelList> newFaceLabelsPtr_;
//- New face labels after mapping
mutable labelList* newFaceLabelsMapPtr_;
mutable std::unique_ptr<labelList> newFaceLabelsMapPtr_;
//- Direct addressing (only one form of addressing is used)
mutable labelList* directAddrPtr_;
mutable std::unique_ptr<labelList> directAddrPtr_;
//- Interpolated addressing (only one form of addressing is used)
mutable labelListList* interpolationAddrPtr_;
mutable std::unique_ptr<labelListList> interpolationAddrPtr_;
//- Interpolation weights
mutable scalarListList* weightsPtr_;
mutable std::unique_ptr<scalarListList> weightsPtr_;
//- Inserted faces
mutable labelList* insertedObjectLabelsPtr_;
mutable std::unique_ptr<labelList> insertedObjectLabelsPtr_;
// Private Member Functions

View File

@ -43,13 +43,13 @@ void Foam::faEdgeMapper::calcAddressing() const
hasUnmapped_ = false;
// Dummy mapping: take value from edge 0
directAddrPtr_ = new labelList(size(), Zero);
directAddrPtr_.reset(new labelList(size(), Zero));
}
void Foam::faEdgeMapper::clearOut()
{
deleteDemandDrivenData(directAddrPtr_);
directAddrPtr_.reset(nullptr);
hasUnmapped_ = false;
}

View File

@ -78,7 +78,7 @@ class faEdgeMapper
mutable bool hasUnmapped_;
//- Direct addressing
mutable labelList* directAddrPtr_;
mutable std::unique_ptr<labelList> directAddrPtr_;
// Private Member Functions

View File

@ -47,7 +47,7 @@ void Foam::faPatchMapper::calcAddressing() const
// Compatibility change HJ, 12/Aug/2017
hasUnmapped_ = false;
directAddrPtr_ = new labelList(patch_.size(), Zero);
directAddrPtr_.reset(new labelList(patch_.size(), Zero));
labelList& addr = *directAddrPtr_;
// Make a map of old edgeFaces, giving edge index in patch given the new
@ -90,7 +90,7 @@ void Foam::faPatchMapper::calcAddressing() const
void Foam::faPatchMapper::clearOut()
{
deleteDemandDrivenData(directAddrPtr_);
directAddrPtr_.reset(nullptr);
hasUnmapped_ = false;
}

View File

@ -84,7 +84,7 @@ class faPatchMapper
mutable bool hasUnmapped_;
//- Direct addressing
mutable labelList* directAddrPtr_;
mutable std::unique_ptr<labelList> directAddrPtr_;
// Private Member Functions

View File

@ -139,7 +139,7 @@ Foam::autoPtr<Foam::faMesh> Foam::faMeshTools::newMesh
IOobject cmptIO(meshIO, "faceLabels", meshSubDir);
cmptIO.readOpt(IOobject::MUST_READ);
cmptIO.writeOpt(IOobject::NO_WRITE);
cmptIO.registerObject(false);
cmptIO.registerObject(IOobject::NO_REGISTER);
// Check who has a mesh

View File

@ -172,11 +172,11 @@ void Foam::faMesh::mapFields(const faMeshMapper& mapper) const
void Foam::faMesh::mapOldAreas(const faMeshMapper& mapper) const
{
if (S0Ptr_)
if (S0ptr_)
{
DebugInFunction << "Mapping old face areas." << endl;
scalarField& S0 = *S0Ptr_;
scalarField& S0 = *S0ptr_;
scalarField savedS0(S0);
S0.setSize(nFaces());
@ -197,11 +197,11 @@ void Foam::faMesh::mapOldAreas(const faMeshMapper& mapper) const
}
}
if (S00Ptr_)
if (S00ptr_)
{
DebugInFunction << "Mapping old-old face areas." << endl;
scalarField& S00 = *S00Ptr_;
scalarField& S00 = *S00ptr_;
scalarField savedS00(S00);
S00.setSize(nFaces());

View File

@ -85,9 +85,9 @@ Foam::wordList Foam::faPatch::constraintTypes()
void Foam::faPatch::clearOut()
{
deleteDemandDrivenData(edgeFacesPtr_);
deleteDemandDrivenData(pointLabelsPtr_);
deleteDemandDrivenData(pointEdgesPtr_);
edgeFacesPtr_.reset(nullptr);
pointLabelsPtr_.reset(nullptr);
pointEdgesPtr_.reset(nullptr);
}
@ -321,7 +321,7 @@ void Foam::faPatch::calcPointLabels() const
}
// Transfer to plain list (reuse storage)
pointLabelsPtr_ = new labelList(std::move(dynEdgePoints));
pointLabelsPtr_.reset(new labelList(std::move(dynEdgePoints)));
/// const labelList& edgePoints = *pointLabelsPtr_;
///
/// // Cannot use invertManyToMany - we have non-local edge numbering
@ -370,7 +370,7 @@ void Foam::faPatch::calcPointEdges() const
}
// Flatten to regular list
pointEdgesPtr_ = new labelListList(edgePoints.size());
pointEdgesPtr_.reset(new labelListList(edgePoints.size()));
auto& pEdges = *pointEdgesPtr_;
forAll(pEdges, pointi)
@ -430,9 +430,12 @@ const Foam::labelUList& Foam::faPatch::edgeFaces() const
{
if (!edgeFacesPtr_)
{
edgeFacesPtr_ = new labelList::subList
edgeFacesPtr_.reset
(
patchSlice(boundaryMesh().mesh().edgeOwner())
new labelList::subList
(
patchSlice(boundaryMesh().mesh().edgeOwner())
)
);
}
@ -478,7 +481,19 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::delta() const
{
// Use patch-normal delta for all non-coupled BCs
const vectorField nHat(edgeNormals());
return nHat*(nHat & (edgeCentres() - edgeFaceCentres()));
vectorField edgePN(edgeCentres() - edgeFaceCentres());
// Do not allow any mag(val) < SMALL
for (vector& edgei : edgePN)
{
if (mag(edgei) < SMALL)
{
edgei = vector::uniform(SMALL);
}
}
return nHat*(nHat & edgePN);
}
@ -490,11 +505,9 @@ void Foam::faPatch::makeDeltaCoeffs(scalarField& dc) const
void Foam::faPatch::makeCorrectionVectors(vectorField& k) const
{
vectorField unitDelta(delta()/mag(delta()));
vectorField edgeNormMag(edgeNormals()/mag(edgeNormals()));
scalarField dn(edgeNormals() & delta());
const vectorField deltaHat(delta()/mag(delta()));
k = edgeNormMag - (scalar(1)/(unitDelta & edgeNormMag))*unitDelta;
k = edgeNormals() - (scalar(1)/(deltaHat & edgeNormals()))*deltaHat;
}

View File

@ -86,13 +86,13 @@ class faPatch
const faBoundaryMesh& boundaryMesh_;
//- Demand-driven: edge-face addressing
mutable labelList::subList* edgeFacesPtr_;
mutable unique_ptr<labelList::subList> edgeFacesPtr_;
//- Demand-driven: local points labels
mutable labelList* pointLabelsPtr_;
mutable unique_ptr<labelList> pointLabelsPtr_;
//- Demand-driven: point-edge addressing
mutable labelListList* pointEdgesPtr_;
mutable unique_ptr<labelListList> pointEdgesPtr_;
// Private Member Functions

View File

@ -84,7 +84,7 @@ void Foam::processorFaPatchField<Foam::scalar>::updateInterfaceMatrix
{
forAll(edgeFaces, facei)
{
result[edgeFaces[facei]] -= coeffs[facei]*pnf[facei];
result[edgeFaces[facei]] += coeffs[facei]*pnf[facei];
}
}
else

View File

@ -48,6 +48,7 @@ SourceFiles
namespace Foam
{
// Forward Declarations
template<class Type>
class faMatrix;
@ -67,8 +68,9 @@ class convectionScheme
:
public refCount
{
// Private data
// Private Data
//- Reference to the mesh
const faMesh& mesh_;
@ -129,11 +131,8 @@ public:
// Member Functions
//- Return mesh reference
const faMesh& mesh() const
{
return mesh_;
}
//- Return const reference to mesh
const faMesh& mesh() const noexcept { return mesh_; }
virtual tmp<GeometricField<Type, faePatchField, edgeMesh>> flux
(

View File

@ -64,13 +64,10 @@ gaussConvectionScheme<Type>::famDiv
tmp<edgeScalarField> tweights = tinterpScheme_().weights(vf);
const edgeScalarField& weights = tweights();
tmp<faMatrix<Type>> tfam
auto tfam = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
faceFlux.dimensions()*vf.dimensions()
)
vf,
faceFlux.dimensions()*vf.dimensions()
);
faMatrix<Type>& fam = tfam.ref();

View File

@ -60,8 +60,9 @@ class gaussConvectionScheme
:
public fa::convectionScheme<Type>
{
// Private data
// Private Data
//- Edge-interpolation scheme
tmp<edgeInterpolationScheme<Type>> tinterpScheme_;

View File

@ -62,40 +62,32 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
const dimensioned<Type> dt
)
{
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_();
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_();
dimensionedScalar rDeltaT2 =
const dimensionedScalar rDeltaT2 =
4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
scalar coefft = (deltaT + deltaT0)/(2*deltaT);
scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
IOobject d2dt2IOobject
const IOobject d2dt2IOobject
(
"d2dt2("+dt.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("d2dt2("+dt.name()+')')
);
if (mesh().moving())
{
scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
const scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
scalarField SS0 = mesh().S() + mesh().S0();
scalarField S0S00 = mesh().S0() + mesh().S00();
const scalarField SS0(mesh().S() + mesh().S0());
const scalarField S0S00(mesh().S0() + mesh().S00());
tmp<GeometricField<Type, faPatchField, areaMesh>> tdt2dt2
auto tdt2dt2 = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
d2dt2IOobject,
mesh(),
dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero)
)
d2dt2IOobject,
mesh(),
dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero)
);
tdt2dt2.ref().primitiveFieldRef() =
@ -107,15 +99,12 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
d2dt2IOobject,
mesh(),
dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero),
calculatedFaPatchField<Type>::typeName
)
d2dt2IOobject,
mesh(),
dimensioned<Type>(dt.dimensions()/dimTime/dimTime, Zero),
calculatedFaPatchField<Type>::typeName
);
}
}
@ -128,70 +117,58 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT2 =
const dimensionedScalar rDeltaT2 =
4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
IOobject d2dt2IOobject
const IOobject d2dt2IOobject
(
"d2dt2("+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("d2dt2("+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_();
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_();
scalar coefft = (deltaT + deltaT0)/(2*deltaT);
scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
scalar coefft0 = coefft + coefft00;
const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
const scalar coefft0 = coefft + coefft00;
if (mesh().moving())
{
scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
const scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
scalarField SS0 = mesh().S() + mesh().S0();
scalarField S0S00 = mesh().S0() + mesh().S00();
const scalarField SS0(mesh().S() + mesh().S0());
const scalarField S0S00(mesh().S0() + mesh().S00());
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
d2dt2IOobject,
mesh(),
rDeltaT2.dimensions()*vf.dimensions(),
halfRdeltaT2*
(
d2dt2IOobject,
mesh(),
rDeltaT2.dimensions()*vf.dimensions(),
halfRdeltaT2*
(
coefft*SS0*vf.primitiveField()
- (coefft*SS0 + coefft00*S0S00)
*vf.oldTime().primitiveField()
+ (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
)/mesh().S(),
rDeltaT2.value()*
(
coefft*vf.boundaryField()
- coefft0*vf.oldTime().boundaryField()
+ coefft00*vf.oldTime().oldTime().boundaryField()
)
coefft*SS0*vf.primitiveField()
- (coefft*SS0 + coefft00*S0S00)
*vf.oldTime().primitiveField()
+ (coefft00*S0S00)*vf.oldTime().oldTime().primitiveField()
)/mesh().S(),
rDeltaT2.value()*
(
coefft*vf.boundaryField()
- coefft0*vf.oldTime().boundaryField()
+ coefft00*vf.oldTime().oldTime().boundaryField()
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
d2dt2IOobject,
rDeltaT2*
(
d2dt2IOobject,
rDeltaT2*
(
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
);
}
@ -206,77 +183,64 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT2 =
const dimensionedScalar rDeltaT2 =
4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
IOobject d2dt2IOobject
const IOobject d2dt2IOobject
(
"d2dt2("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("d2dt2("+rho.name()+','+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_();
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_();
scalar coefft = (deltaT + deltaT0)/(2*deltaT);
scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
scalar coefft0 = coefft + coefft00;
const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
const scalar coefft0 = coefft + coefft00;
if (mesh().moving())
{
scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
const scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
scalarField SS0rhoRho0 =
(mesh().S() + mesh().S0())
*rho.value();
const scalarField SS0rhoRho0((mesh().S() + mesh().S0())*rho.value());
scalarField S0S00rho0Rho00 =
(mesh().S0() + mesh().S00())
*rho.value();
return tmp<GeometricField<Type, faPatchField, areaMesh>>
const scalarField S0S00rho0Rho00
(
new GeometricField<Type, faPatchField, areaMesh>
(mesh().S0() + mesh().S00())
*rho.value()
);
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
d2dt2IOobject,
mesh(),
rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
halfRdeltaT2*
(
d2dt2IOobject,
mesh(),
rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
halfRdeltaT2*
(
coefft*SS0rhoRho0*vf.primitiveField()
- (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
*vf.oldTime().primitiveField()
+ (coefft00*S0S00rho0Rho00)
*vf.oldTime().oldTime().primitiveField()
)/mesh().S(),
rDeltaT2.value()*rho.value()*
(
coefft*vf.boundaryField()
- coefft0*vf.oldTime().boundaryField()
+ coefft00*vf.oldTime().oldTime().boundaryField()
)
coefft*SS0rhoRho0*vf.primitiveField()
- (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
*vf.oldTime().primitiveField()
+ (coefft00*S0S00rho0Rho00)
*vf.oldTime().oldTime().primitiveField()
)/mesh().S(),
rDeltaT2.value()*rho.value()*
(
coefft*vf.boundaryField()
- coefft0*vf.oldTime().boundaryField()
+ coefft00*vf.oldTime().oldTime().boundaryField()
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
d2dt2IOobject,
rDeltaT2 * rho.value() *
(
d2dt2IOobject,
rDeltaT2 * rho.value() *
(
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
);
}
@ -291,36 +255,32 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT2 =
const dimensionedScalar rDeltaT2 =
4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());
IOobject d2dt2IOobject
const IOobject d2dt2IOobject
(
"d2dt2("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("d2dt2("+rho.name()+','+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_();
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_();
scalar coefft = (deltaT + deltaT0)/(2*deltaT);
scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
if (mesh().moving())
{
scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
const scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
const scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();
scalarField SS0rhoRho0
const scalarField SS0rhoRho0
(
(mesh().S() + mesh().S0())
*(rho.primitiveField() + rho.oldTime().primitiveField())
);
scalarField S0S00rho0Rho00
const scalarField S0S00rho0Rho00
(
(mesh().S0() + mesh().S00())
*(
@ -329,69 +289,63 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
)
);
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
d2dt2IOobject,
mesh(),
rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
quarterRdeltaT2*
(
d2dt2IOobject,
mesh(),
rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
quarterRdeltaT2*
(
coefft*SS0rhoRho0*vf.primitiveField()
coefft*SS0rhoRho0*vf.primitiveField()
- (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
*vf.oldTime().primitiveField()
- (coefft*SS0rhoRho0 + coefft00*S0S00rho0Rho00)
*vf.oldTime().primitiveField()
+ (coefft00*S0S00rho0Rho00)
*vf.oldTime().oldTime().primitiveField()
)/mesh().S(),
halfRdeltaT2*
(
+ (coefft00*S0S00rho0Rho00)
*vf.oldTime().oldTime().primitiveField()
)/mesh().S(),
halfRdeltaT2*
(
coefft
*(rho.boundaryField() + rho.oldTime().boundaryField())
*vf.boundaryField()
- (
coefft
*(rho.boundaryField() + rho.oldTime().boundaryField())
*vf.boundaryField()
- (
coefft
*(
rho.boundaryField()
+ rho.oldTime().boundaryField()
)
+ coefft00
*(
rho.oldTime().boundaryField()
+ rho.oldTime().oldTime().boundaryField()
)
)*vf.oldTime().boundaryField()
*(
rho.boundaryField()
+ rho.oldTime().boundaryField()
)
+ coefft00
*(
rho.oldTime().boundaryField()
+ rho.oldTime().oldTime().boundaryField()
)*vf.oldTime().oldTime().boundaryField()
)
rho.oldTime().boundaryField()
+ rho.oldTime().oldTime().boundaryField()
)
)*vf.oldTime().boundaryField()
+ coefft00
*(
rho.oldTime().boundaryField()
+ rho.oldTime().oldTime().boundaryField()
)*vf.oldTime().oldTime().boundaryField()
)
);
}
else
{
dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
const dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;
areaScalarField rhoRho0(rho + rho.oldTime());
areaScalarField rho0Rho00(rho.oldTime() + rho.oldTime().oldTime());
const areaScalarField rhoRho0(rho + rho.oldTime());
const areaScalarField rho0Rho00(rho.oldTime() + rho.oldTime().oldTime());
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
d2dt2IOobject,
halfRdeltaT2*
(
d2dt2IOobject,
halfRdeltaT2*
(
coefft*rhoRho0*vf
- (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
+ coefft00*rho0Rho00*vf.oldTime().oldTime()
)
coefft*rhoRho0*vf
- (coefft*rhoRho0 + coefft00*rho0Rho00)*vf.oldTime()
+ coefft00*rho0Rho00*vf.oldTime().oldTime()
)
);
}
@ -405,32 +359,28 @@ EulerFaD2dt2Scheme<Type>::famD2dt2
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
auto tfam = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
vf.dimensions()*dimArea/dimTime/dimTime
)
vf,
vf.dimensions()*dimArea/dimTime/dimTime
);
faMatrix<Type>& fam = tfam.ref();
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_();
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_();
scalar coefft = (deltaT + deltaT0)/(2*deltaT);
scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
scalar coefft0 = coefft + coefft00;
const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
const scalar coefft0 = coefft + coefft00;
scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
const scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
if (mesh().moving())
{
scalar halfRdeltaT2 = rDeltaT2/2.0;
const scalar halfRdeltaT2 = rDeltaT2/2.0;
scalarField SS0(mesh().S() + mesh().S0());
scalarField S0S00(mesh().S0() + mesh().S00());
const scalarField SS0(mesh().S() + mesh().S0());
const scalarField S0S00(mesh().S0() + mesh().S00());
fam.diag() = (coefft*halfRdeltaT2)*SS0;
@ -465,33 +415,27 @@ EulerFaD2dt2Scheme<Type>::famD2dt2
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
auto tfam = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
rho.dimensions()*vf.dimensions()*dimArea
/dimTime/dimTime
)
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime/dimTime
);
faMatrix<Type>& fam = tfam.ref();
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_();
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_();
scalar coefft = (deltaT + deltaT0)/(2*deltaT);
scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
const scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
if (mesh().moving())
{
scalar halfRdeltaT2 = 0.5*rDeltaT2;
const scalar halfRdeltaT2 = 0.5*rDeltaT2;
scalarField SS0(mesh().S() + mesh().S0());
scalarField S0S00(mesh().S0() + mesh().S00());
const scalarField SS0(mesh().S() + mesh().S0());
const scalarField S0S00(mesh().S0() + mesh().S00());
fam.diag() = rho.value()*(coefft*halfRdeltaT2)*SS0;
@ -526,35 +470,32 @@ EulerFaD2dt2Scheme<Type>::famD2dt2
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
auto tfam = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime/dimTime
)
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime/dimTime
);
faMatrix<Type>& fam = tfam.ref();
scalar deltaT =deltaT_();
scalar deltaT0 = deltaT0_();
const scalar deltaT =deltaT_();
const scalar deltaT0 = deltaT0_();
scalar coefft = (deltaT + deltaT0)/(2*deltaT);
scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
const scalar coefft = (deltaT + deltaT0)/(2*deltaT);
const scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);
scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
const scalar rDeltaT2 = 4.0/sqr(deltaT + deltaT0);
if (mesh().moving())
{
scalar quarterRdeltaT2 = 0.25*rDeltaT2;
const scalar quarterRdeltaT2 = 0.25*rDeltaT2;
scalarField SS0rhoRho0
const scalarField SS0rhoRho0
(
(mesh().S() + mesh().S0())
*(rho.primitiveField() + rho.oldTime().primitiveField())
);
scalarField S0S00rho0Rho00
const scalarField S0S00rho0Rho00
(
(mesh().S0() + mesh().S00())
*(
@ -576,14 +517,14 @@ EulerFaD2dt2Scheme<Type>::famD2dt2
}
else
{
scalar halfRdeltaT2 = 0.5*rDeltaT2;
const scalar halfRdeltaT2 = 0.5*rDeltaT2;
scalarField rhoRho0
const scalarField rhoRho0
(
rho.primitiveField() + rho.oldTime().primitiveField()
);
scalarField rho0Rho00
const scalarField rho0Rho00
(
rho.oldTime().primitiveField()
+ rho.oldTime().oldTime().primitiveField()

View File

@ -98,7 +98,7 @@ public:
// Member Functions
//- Return mesh reference
const faMesh& mesh() const
const faMesh& mesh() const noexcept
{
return fa::faD2dt2Scheme<Type>::mesh();
}

View File

@ -24,13 +24,11 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Abstract base class for finite area d2dt2 schemes.
\*---------------------------------------------------------------------------*/
#include "fa.H"
#include "HashTable.H"
#include "IOobject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -42,6 +40,24 @@ namespace Foam
namespace fa
{
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class Type>
IOobject faD2dt2Scheme<Type>::fieldIOobject(const word& fldName) const
{
const auto& pMesh = mesh().mesh();
return IOobject
(
fldName,
pMesh.time().timeName(),
pMesh,
IOobject::NO_READ,
IOobject::NO_WRITE
);
}
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
template<class Type>

View File

@ -49,10 +49,12 @@ SourceFiles
namespace Foam
{
// Forward Declarations
template<class Type>
class faMatrix;
class faMesh;
class IOobject;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -68,15 +70,18 @@ class faD2dt2Scheme
:
public refCount
{
protected:
// Protected data
// Protected Data
//- Reference to the mesh
const faMesh& mesh_;
// Private Member Functions
// Protected Member Functions
//- Helper to construct IOobject for field and current time
IOobject fieldIOobject(const word& fldName) const;
//- No copy construct
faD2dt2Scheme(const faD2dt2Scheme&) = delete;
@ -135,10 +140,7 @@ public:
// Member Functions
//- Return mesh reference
const faMesh& mesh() const
{
return mesh_;
}
const faMesh& mesh() const noexcept { return mesh_; }
virtual tmp<GeometricField<Type, faPatchField, areaMesh>> facD2dt2
(

View File

@ -48,27 +48,20 @@ EulerFaDdtScheme<Type>::facDdt
const dimensioned<Type> dt
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+dt.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+dt.name()+')')
);
if (mesh().moving())
{
tmp<GeometricField<Type, faPatchField, areaMesh>> tdtdt
auto tdtdt = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
mesh(),
dimensioned<Type>(dt.dimensions()/dimTime, Zero)
)
ddtIOobject,
mesh(),
dimensioned<Type>(dt.dimensions()/dimTime, Zero)
);
tdtdt.ref().primitiveFieldRef() =
@ -95,25 +88,18 @@ EulerFaDdtScheme<Type>::facDdt0
const dimensioned<Type> dt
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+dt.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+dt.name()+')')
);
tmp<GeometricField<Type, faPatchField, areaMesh>> tdtdt0
auto tdtdt0 = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
mesh(),
-rDeltaT*dt
)
ddtIOobject,
mesh(),
-rDeltaT*dt
);
if (mesh().moving())
@ -133,47 +119,37 @@ EulerFaDdtScheme<Type>::facDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+vf.name()+')')
);
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
vf()
- vf.oldTime()()*mesh().S0()/mesh().S()
),
rDeltaT.value()*
(
vf.boundaryField() - vf.oldTime().boundaryField()
)
vf()
- vf.oldTime()()*mesh().S0()/mesh().S()
),
rDeltaT.value()*
(
vf.boundaryField() - vf.oldTime().boundaryField()
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
rDeltaT*(vf - vf.oldTime())
)
ddtIOobject,
rDeltaT*(vf - vf.oldTime())
);
}
}
@ -186,40 +162,30 @@ EulerFaDdtScheme<Type>::facDdt0
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt0("+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt0("+vf.name()+')')
);
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
(-rDeltaT.value())*vf.oldTime().internalField(),
(-rDeltaT.value())*vf.oldTime().boundaryField()
)
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
(-rDeltaT.value())*vf.oldTime().internalField(),
(-rDeltaT.value())*vf.oldTime().boundaryField()
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
(-rDeltaT)*vf.oldTime()
)
ddtIOobject,
(-rDeltaT)*vf.oldTime()
);
}
}
@ -233,47 +199,37 @@ EulerFaDdtScheme<Type>::facDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')')
);
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*rho.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*rho.value()*
(
vf()
- vf.oldTime()()*mesh().S0()/mesh().S()
),
rDeltaT.value()*rho.value()*
(
vf.boundaryField() - vf.oldTime().boundaryField()
)
vf()
- vf.oldTime()()*mesh().S0()/mesh().S()
),
rDeltaT.value()*rho.value()*
(
vf.boundaryField() - vf.oldTime().boundaryField()
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
rDeltaT*rho*(vf - vf.oldTime())
)
ddtIOobject,
rDeltaT*rho*(vf - vf.oldTime())
);
}
}
@ -286,42 +242,32 @@ EulerFaDdtScheme<Type>::facDdt0
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt0("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')')
);
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
(-rDeltaT.value())*rho.value()*
vf.oldTime()()*mesh().S0()/mesh().S(),
(-rDeltaT.value())*rho.value()*
vf.oldTime().boundaryField()
)
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
(-rDeltaT.value())*rho.value()*
vf.oldTime()()*mesh().S0()/mesh().S(),
(-rDeltaT.value())*rho.value()*
vf.oldTime().boundaryField()
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
(-rDeltaT)*rho*vf.oldTime()
)
ddtIOobject,
(-rDeltaT)*rho*vf.oldTime()
);
}
}
@ -335,50 +281,40 @@ EulerFaDdtScheme<Type>::facDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')')
);
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
rho()*vf()
- rho.oldTime()()
*vf.oldTime()()*mesh().S0()/mesh().S()
),
rDeltaT.value()*
(
rho.boundaryField()*vf.boundaryField()
- rho.oldTime().boundaryField()
*vf.oldTime().boundaryField()
)
rho()*vf()
- rho.oldTime()()
*vf.oldTime()()*mesh().S0()/mesh().S()
),
rDeltaT.value()*
(
rho.boundaryField()*vf.boundaryField()
- rho.oldTime().boundaryField()
*vf.oldTime().boundaryField()
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
)
ddtIOobject,
rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
);
}
}
@ -392,52 +328,43 @@ EulerFaDdtScheme<Type>::facDdt0
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt0("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')')
);
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
- rho.oldTime()()
*vf.oldTime()()*mesh().S0()/mesh().S()
),
rDeltaT.value()*
(
- rho.oldTime().boundaryField()
*vf.oldTime().boundaryField()
)
- rho.oldTime()()
*vf.oldTime()()*mesh().S0()/mesh().S()
),
rDeltaT.value()*
(
- rho.oldTime().boundaryField()
*vf.oldTime().boundaryField()
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
(-rDeltaT)*rho.oldTime()*vf.oldTime()
)
ddtIOobject,
(-rDeltaT)*rho.oldTime()*vf.oldTime()
);
}
}
template<class Type>
tmp<faMatrix<Type>>
EulerFaDdtScheme<Type>::famDdt
@ -445,18 +372,14 @@ EulerFaDdtScheme<Type>::famDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
auto tfam = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
vf.dimensions()*dimArea/dimTime
)
vf,
vf.dimensions()*dimArea/dimTime
);
faMatrix<Type>& fam = tfam.ref();
scalar rDeltaT = 1.0/mesh().time().deltaT().value();
const scalar rDeltaT = 1.0/mesh().time().deltaT().value();
fam.diag() = rDeltaT*mesh().S();
@ -481,17 +404,14 @@ EulerFaDdtScheme<Type>::famDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
auto tfam = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
)
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
);
faMatrix<Type>& fam = tfam.ref();
scalar rDeltaT = 1.0/mesh().time().deltaT().value();
const scalar rDeltaT = 1.0/mesh().time().deltaT().value();
fam.diag() = rDeltaT*rho.value()*mesh().S();
@ -518,17 +438,14 @@ EulerFaDdtScheme<Type>::famDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
auto tfam = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
)
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
);
faMatrix<Type>& fam = tfam.ref();
scalar rDeltaT = 1.0/mesh().time().deltaT().value();
const scalar rDeltaT = 1.0/mesh().time().deltaT().value();
fam.diag() = rDeltaT*rho.primitiveField()*mesh().S();

View File

@ -63,10 +63,8 @@ scalar backwardFaDdtScheme<Type>::deltaT0_(const GeoField& vf) const
{
return GREAT;
}
else
{
return deltaT0_();
}
return deltaT0_();
}
@ -79,34 +77,27 @@ backwardFaDdtScheme<Type>::facDdt
const dimensioned<Type> dt
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+dt.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+dt.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_();
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_();
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
if (mesh().moving())
{
tmp<GeometricField<Type, faPatchField, areaMesh>> tdtdt
auto tdtdt = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
mesh(),
dimensioned<Type>(dt.dimensions()/dimTime, Zero)
)
ddtIOobject,
mesh(),
dimensioned<Type>(dt.dimensions()/dimTime, Zero)
);
tdtdt.ref().primitiveFieldRef() = rDeltaT.value()*dt.value()*
@ -134,32 +125,25 @@ backwardFaDdtScheme<Type>::facDdt0
const dimensioned<Type> dt
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+dt.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+dt.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_();
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_();
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
tmp<GeometricField<Type, faPatchField, areaMesh>> tdtdt0
auto tdtdt0 = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
(
ddtIOobject,
mesh(),
-rDeltaT*(coefft0 - coefft00)*dt
)
ddtIOobject,
mesh(),
-rDeltaT*(coefft0 - coefft00)*dt
);
if (mesh().moving())
@ -181,66 +165,56 @@ backwardFaDdtScheme<Type>::facDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
rDeltaT.value()*
coefft*vf() -
(
coefft*vf() -
(
coefft0*vf.oldTime()()*mesh().S0()
- coefft00*vf.oldTime().oldTime()()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
coefft0*vf.oldTime()()*mesh().S0()
- coefft00*vf.oldTime().oldTime()()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
coefft*vf.boundaryField() -
(
coefft*vf.boundaryField() -
(
coefft0*vf.oldTime().boundaryField()
- coefft00*vf.oldTime().oldTime().boundaryField()
)
coefft0*vf.oldTime().boundaryField()
- coefft00*vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
rDeltaT*
(
ddtIOobject,
rDeltaT*
(
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
);
}
@ -254,63 +228,53 @@ backwardFaDdtScheme<Type>::facDdt0
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt0("+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt0("+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
- (
coefft0*vf.oldTime()()*mesh().S0()
- coefft00*vf.oldTime().oldTime()()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
- (
coefft0*vf.oldTime().boundaryField()
- coefft00*vf.oldTime().oldTime().boundaryField()
)
- (
coefft0*vf.oldTime()()*mesh().S0()
- coefft00*vf.oldTime().oldTime()()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
- (
coefft0*vf.oldTime().boundaryField()
- coefft00*vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
rDeltaT*
(
ddtIOobject,
rDeltaT*
(
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
);
}
@ -325,71 +289,62 @@ backwardFaDdtScheme<Type>::facDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*rho.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*rho.value()*
coefft*vf.internalField() -
(
coefft*vf.internalField() -
(
coefft0*vf.oldTime()()*mesh().S0()
- coefft00*vf.oldTime().oldTime()()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*rho.value()*
coefft0*vf.oldTime()()*mesh().S0()
- coefft00*vf.oldTime().oldTime()()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*rho.value()*
(
coefft*vf.boundaryField() -
(
coefft*vf.boundaryField() -
(
coefft0*vf.oldTime().boundaryField()
- coefft00*vf.oldTime().oldTime().boundaryField()
)
coefft0*vf.oldTime().boundaryField()
- coefft00*vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
rDeltaT*rho*
(
ddtIOobject,
rDeltaT*rho*
(
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
);
}
}
template<class Type>
tmp<GeometricField<Type, faPatchField, areaMesh>>
backwardFaDdtScheme<Type>::facDdt0
@ -398,63 +353,53 @@ backwardFaDdtScheme<Type>::facDdt0
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt0("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*rho.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*rho.value()*
(
-(
coefft0*vf.oldTime()()*mesh().S0()
- coefft00*vf.oldTime().oldTime()()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*rho.value()*
(
-(
coefft0*vf.oldTime().boundaryField()
- coefft00*vf.oldTime().oldTime().boundaryField()
)
- (
coefft0*vf.oldTime()()*mesh().S0()
- coefft00*vf.oldTime().oldTime()()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*rho.value()*
(
- (
coefft0*vf.oldTime().boundaryField()
- coefft00*vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
rDeltaT*rho*
(
ddtIOobject,
rDeltaT*rho*
(
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
);
}
@ -469,69 +414,59 @@ backwardFaDdtScheme<Type>::facDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
coefft*rho.internalField()*vf.internalField() -
(
coefft*rho.internalField()*vf.internalField() -
(
coefft0*rho.oldTime()()
*vf.oldTime()()*mesh().S0()
- coefft00*rho.oldTime().oldTime()()
*vf.oldTime().oldTime()()*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
coefft0*rho.oldTime()()
*vf.oldTime()()*mesh().S0()
- coefft00*rho.oldTime().oldTime()()
*vf.oldTime().oldTime()()*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
coefft*vf.boundaryField() -
(
coefft*vf.boundaryField() -
(
coefft0*rho.oldTime().boundaryField()
*vf.oldTime().boundaryField()
- coefft00*rho.oldTime().oldTime().boundaryField()
*vf.oldTime().oldTime().boundaryField()
)
coefft0*rho.oldTime().boundaryField()
*vf.oldTime().boundaryField()
- coefft00*rho.oldTime().oldTime().boundaryField()
*vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
rDeltaT*
(
ddtIOobject,
rDeltaT*
(
coefft*rho*vf
- coefft0*rho.oldTime()*vf.oldTime()
+ coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
)
coefft*rho*vf
- coefft0*rho.oldTime()*vf.oldTime()
+ coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
)
);
}
@ -546,66 +481,56 @@ backwardFaDdtScheme<Type>::facDdt0
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt0("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
if (mesh().moving())
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
- (
coefft0*rho.oldTime()()
*vf.oldTime()()*mesh().S0()
- coefft00*rho.oldTime().oldTime()()
*vf.oldTime().oldTime()()*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
- (
coefft0*rho.oldTime().boundaryField()
*vf.oldTime().boundaryField()
- coefft00*rho.oldTime().oldTime().boundaryField()
*vf.oldTime().oldTime().boundaryField()
)
- (
coefft0*rho.oldTime()()
*vf.oldTime()()*mesh().S0()
- coefft00*rho.oldTime().oldTime()()
*vf.oldTime().oldTime()()*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
- (
coefft0*rho.oldTime().boundaryField()
*vf.oldTime().boundaryField()
- coefft00*rho.oldTime().oldTime().boundaryField()
*vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
ddtIOobject,
rDeltaT*
(
ddtIOobject,
rDeltaT*
(
- coefft0*rho.oldTime()*vf.oldTime()
+ coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
)
- coefft0*rho.oldTime()*vf.oldTime()
+ coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
)
);
}
@ -619,25 +544,22 @@ backwardFaDdtScheme<Type>::famDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
auto tfam = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
vf.dimensions()*dimArea/dimTime
)
vf,
vf.dimensions()*dimArea/dimTime
);
faMatrix<Type>& fam = tfam.ref();
scalar rDeltaT = 1.0/deltaT_();
const scalar rDeltaT = 1.0/deltaT_();
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
fam.diag() = (coefft*rDeltaT)*mesh().S();
@ -671,24 +593,21 @@ backwardFaDdtScheme<Type>::famDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
auto tfam = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
)
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
);
faMatrix<Type>& fam = tfam.ref();
scalar rDeltaT = 1.0/deltaT_();
const scalar rDeltaT = 1.0/deltaT_();
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
fam.diag() = (coefft*rDeltaT*rho.value())*mesh().S();
@ -722,24 +641,21 @@ backwardFaDdtScheme<Type>::famDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
auto tfam = tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
)
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
);
faMatrix<Type>& fam = tfam.ref();
scalar rDeltaT = 1.0/deltaT_();
const scalar rDeltaT = 1.0/deltaT_();
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
fam.diag() = (coefft*rDeltaT)*rho.primitiveField()*mesh().S();

View File

@ -62,34 +62,27 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
{
// No change compared to backward
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+dt.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+dt.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_();
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_();
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
if (mesh().moving())
{
tmp<areaScalarField> tdtdt
auto tdtdt = tmp<areaScalarField>::New
(
new areaScalarField
(
ddtIOobject,
mesh(),
dimensionedScalar(dt.dimensions()/dimTime, Zero)
)
ddtIOobject,
mesh(),
dimensionedScalar(dt.dimensions()/dimTime, Zero)
);
tdtdt.ref().primitiveFieldRef() = rDeltaT.value()*dt.value()*
@ -117,32 +110,25 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
{
// No change compared to backward
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+dt.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+dt.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_();
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_();
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
const scalar coefft = 1 + deltaT/(deltaT + deltaT0);
const scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
const scalar coefft0 = coefft + coefft00;
tmp<areaScalarField> tdtdt0
auto tdtdt0 = tmp<areaScalarField>::New
(
new areaScalarField
(
ddtIOobject,
mesh(),
-rDeltaT*(coefft0 - coefft00)*dt
)
ddtIOobject,
mesh(),
-rDeltaT*(coefft0 - coefft00)*dt
);
if (mesh().moving())
@ -162,24 +148,20 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
const areaScalarField& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
// Calculate unboundedness indicator
// Note: all times moved by one because access to internal field
// copies current field into the old-time level.
areaScalarField phict
const areaScalarField phict
(
mag
(
@ -192,62 +174,59 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
vf.oldTime()
- vf.oldTime().oldTime()
)
+ dimensionedScalar("small", vf.dimensions(), SMALL)
+ dimensionedScalar(vf.dimensions(), SMALL)
)
);
areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
areaScalarField coefft0(coefft + coefft00);
const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
const areaScalarField coefft00
(
limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
);
const areaScalarField coefft0(coefft + coefft00);
if (mesh().moving())
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
rDeltaT.value()*
coefft*vf.primitiveField() -
(
coefft*vf.primitiveField() -
(
coefft0.primitiveField()
*vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()
*vf.oldTime().oldTime().primitiveField()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
coefft0.primitiveField()
*vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()
*vf.oldTime().oldTime().primitiveField()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
coefft.boundaryField()*vf.boundaryField() -
(
coefft.boundaryField()*vf.boundaryField() -
(
coefft0.boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
coefft0.boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
rDeltaT*
(
ddtIOobject,
rDeltaT*
(
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
);
}
@ -259,24 +238,20 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
const areaScalarField& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt0("+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt0("+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
// Calculate unboundedness indicator
// Note: all times moved by one because access to internal field
// copies current field into the old-time level.
areaScalarField phict
const areaScalarField phict
(
mag
(
@ -289,59 +264,56 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
vf.oldTime()
- vf.oldTime().oldTime()
)
+ dimensionedScalar("small", vf.dimensions(), SMALL)
+ dimensionedScalar(vf.dimensions(), SMALL)
)
);
areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
areaScalarField coefft0(coefft + coefft00);
const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
const areaScalarField coefft00
(
limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
);
const areaScalarField coefft0(coefft + coefft00);
if (mesh().moving())
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
- (
coefft0.primitiveField()*
vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()*
vf.oldTime().oldTime().primitiveField()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
- (
coefft0.boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
- (
coefft0.primitiveField()*
vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()*
vf.oldTime().oldTime().primitiveField()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
- (
coefft0.boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
rDeltaT*
(
ddtIOobject,
rDeltaT*
(
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
);
}
@ -354,24 +326,20 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
const areaScalarField& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
// Calculate unboundedness indicator
// Note: all times moved by one because access to internal field
// copies current field into the old-time level.
areaScalarField phict
const areaScalarField phict
(
mag
(
@ -388,87 +356,81 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
)
);
areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
areaScalarField coefft0(coefft + coefft00);
const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
const areaScalarField coefft00
(
limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
);
const areaScalarField coefft0(coefft + coefft00);
if (mesh().moving())
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*rho.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*rho.value()*
coefft*vf.primitiveField() -
(
coefft*vf.primitiveField() -
(
coefft0.primitiveField()*
vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()*
vf.oldTime().oldTime().primitiveField()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*rho.value()*
coefft0.primitiveField()*
vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()*
vf.oldTime().oldTime().primitiveField()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*rho.value()*
(
coefft.boundaryField()*vf.boundaryField() -
(
coefft.boundaryField()*vf.boundaryField() -
(
coefft0.boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
coefft0.boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
rDeltaT*rho*
(
ddtIOobject,
rDeltaT*rho*
(
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
coefft*vf
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
);
}
}
tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
(
const dimensionedScalar& rho,
const areaScalarField& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt0("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
// Calculate unboundedness indicator
// Note: all times moved by one because access to internal field
// copies current field into the old-time level.
areaScalarField phict
const areaScalarField phict
(
mag
(
@ -481,59 +443,56 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
vf.oldTime()
- vf.oldTime().oldTime()
)
+ dimensionedScalar("small", vf.dimensions(), SMALL)
+ dimensionedScalar(vf.dimensions(), SMALL)
)
);
areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
areaScalarField coefft0(coefft + coefft00);
const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
const areaScalarField coefft00
(
limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
);
const areaScalarField coefft0(coefft + coefft00);
if (mesh().moving())
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*rho.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*rho.value()*
(
-(
coefft0.primitiveField()*
vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()*
vf.oldTime().oldTime().primitiveField()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*rho.value()*
(
-(
coefft0.boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
-(
coefft0.primitiveField()*
vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()*
vf.oldTime().oldTime().primitiveField()
*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*rho.value()*
(
-(
coefft0.boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
rDeltaT*rho*
(
ddtIOobject,
rDeltaT*rho*
(
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
- coefft0*vf.oldTime()
+ coefft00*vf.oldTime().oldTime()
)
);
}
@ -546,24 +505,20 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
const areaScalarField& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
// Calculate unboundedness indicator
// Note: all times moved by one because access to internal field
// copies current field into the old-time level.
areaScalarField phict
const areaScalarField phict
(
mag
(
@ -576,65 +531,62 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt
rho.oldTime()*vf.oldTime()
- rho.oldTime().oldTime()*vf.oldTime().oldTime()
)
+ dimensionedScalar("small", rho.dimensions()*vf.dimensions(), SMALL)
+ dimensionedScalar(rho.dimensions()*vf.dimensions(), SMALL)
)
);
areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
areaScalarField coefft0(coefft + coefft00);
const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
const areaScalarField coefft00
(
limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
);
const areaScalarField coefft0(coefft + coefft00);
if (mesh().moving())
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
coefft*rho.primitiveField()*vf.primitiveField() -
(
coefft*rho.primitiveField()*vf.primitiveField() -
(
coefft0.primitiveField()*
rho.oldTime().primitiveField()*
vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()*
rho.oldTime().oldTime().primitiveField()
*vf.oldTime().oldTime().primitiveField()*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
coefft0.primitiveField()*
rho.oldTime().primitiveField()*
vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()*
rho.oldTime().oldTime().primitiveField()
*vf.oldTime().oldTime().primitiveField()*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
coefft.boundaryField()*vf.boundaryField() -
(
coefft.boundaryField()*vf.boundaryField() -
(
coefft0.boundaryField()*
rho.oldTime().boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
rho.oldTime().oldTime().boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
coefft0.boundaryField()*
rho.oldTime().boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
rho.oldTime().oldTime().boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
rDeltaT*
(
ddtIOobject,
rDeltaT*
(
coefft*rho*vf
- coefft0*rho.oldTime()*vf.oldTime()
+ coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
)
coefft*rho*vf
- coefft0*rho.oldTime()*vf.oldTime()
+ coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
)
);
}
@ -647,24 +599,20 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
const areaScalarField& vf
)
{
dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
const dimensionedScalar rDeltaT(1.0/mesh().time().deltaT());
IOobject ddtIOobject
const IOobject ddtIOobject
(
"ddt0("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE
this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')')
);
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
// Calculate unboundedness indicator
// Note: all times moved by one because access to internal field
// copies current field into the old-time level.
areaScalarField phict
const areaScalarField phict
(
mag
(
@ -677,62 +625,59 @@ tmp<areaScalarField> boundedBackwardFaDdtScheme::facDdt0
rho.oldTime()*vf.oldTime()
- rho.oldTime().oldTime()*vf.oldTime().oldTime()
)
+ dimensionedScalar("small", rho.dimensions()*vf.dimensions(), SMALL)
+ dimensionedScalar(rho.dimensions()*vf.dimensions(), SMALL)
)
);
areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
const areaScalarField limiter(pos(phict) - pos(phict - scalar(1)));
areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
areaScalarField coefft00(limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0)));
areaScalarField coefft0(coefft + coefft00);
const areaScalarField coefft(scalar(1) + limiter*deltaT/(deltaT + deltaT0));
const areaScalarField coefft00
(
limiter*sqr(deltaT)/(deltaT0*(deltaT + deltaT0))
);
const areaScalarField coefft0(coefft + coefft00);
if (mesh().moving())
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
ddtIOobject,
mesh(),
rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
rDeltaT.value()*
(
- (
coefft0.primitiveField()*
rho.oldTime().primitiveField()*
vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()*
rho.oldTime().oldTime().primitiveField()*
vf.oldTime().oldTime().primitiveField()*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
- (
coefft0.boundaryField()*
rho.oldTime().boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
rho.oldTime().oldTime().boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
- (
coefft0.primitiveField()*
rho.oldTime().primitiveField()*
vf.oldTime().primitiveField()*mesh().S0()
- coefft00.primitiveField()*
rho.oldTime().oldTime().primitiveField()*
vf.oldTime().oldTime().primitiveField()*mesh().S00()
)/mesh().S()
),
rDeltaT.value()*
(
- (
coefft0.boundaryField()*
rho.oldTime().boundaryField()*
vf.oldTime().boundaryField()
- coefft00.boundaryField()*
rho.oldTime().oldTime().boundaryField()*
vf.oldTime().oldTime().boundaryField()
)
)
);
}
else
{
return tmp<areaScalarField>
return tmp<areaScalarField>::New
(
new areaScalarField
ddtIOobject,
rDeltaT*
(
ddtIOobject,
rDeltaT*
(
- coefft0*rho.oldTime()*vf.oldTime()
+ coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
)
- coefft0*rho.oldTime()*vf.oldTime()
+ coefft00*rho.oldTime().oldTime()*vf.oldTime().oldTime()
)
);
}
@ -744,26 +689,22 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
const areaScalarField& vf
)
{
tmp<faScalarMatrix> tfam
auto tfam = tmp<faScalarMatrix>::New
(
new faScalarMatrix
(
vf,
vf.dimensions()*dimArea/dimTime
)
vf,
vf.dimensions()*dimArea/dimTime
);
faScalarMatrix& fam = tfam.ref();
scalar rDeltaT = 1.0/deltaT_();
const scalar rDeltaT = 1.0/deltaT_();
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
// Calculate unboundedness indicator
// Note: all times moved by one because access to internal field
// copies current field into the old-time level.
scalarField phict
const scalarField phict
(
mag
(
@ -780,11 +721,14 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
)
);
scalarField limiter(pos(phict) - pos(phict - 1.0));
const scalarField limiter(pos(phict) - pos(phict - 1.0));
scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)));
scalarField coefft0(coefft + coefft00);
const scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
const scalarField coefft00
(
limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))
);
const scalarField coefft0(coefft + coefft00);
fam.diag() = (coefft*rDeltaT)*mesh().S();
@ -816,25 +760,22 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
const areaScalarField& vf
)
{
tmp<faScalarMatrix> tfam
auto tfam = tmp<faScalarMatrix>::New
(
new faScalarMatrix
(
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
)
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
);
faScalarMatrix& fam = tfam.ref();
scalar rDeltaT = 1.0/deltaT_();
const scalar rDeltaT = 1.0/deltaT_();
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
// Calculate unboundedness indicator
// Note: all times moved by one because access to internal field
// copies current field into the old-time level.
scalarField phict
const scalarField phict
(
mag
(
@ -851,11 +792,14 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
)
);
scalarField limiter(pos(phict) - pos(phict - 1.0));
const scalarField limiter(pos(phict) - pos(phict - 1.0));
scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)));
scalarField coefft0(coefft + coefft00);
const scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
const scalarField coefft00
(
limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))
);
const scalarField coefft0(coefft + coefft00);
fam.diag() = (coefft*rDeltaT*rho.value())*mesh().S();
@ -887,25 +831,22 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
const areaScalarField& vf
)
{
tmp<faScalarMatrix> tfam
auto tfam = tmp<faScalarMatrix>::New
(
new faScalarMatrix
(
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
)
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
);
faScalarMatrix& fam = tfam.ref();
scalar rDeltaT = 1.0/deltaT_();
const scalar rDeltaT = 1.0/deltaT_();
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
const scalar deltaT = deltaT_();
const scalar deltaT0 = deltaT0_(vf);
// Calculate unboundedness indicator
// Note: all times moved by one because access to internal field
// copies current field into the old-time level.
scalarField phict
const scalarField phict
(
mag
(
@ -926,11 +867,14 @@ tmp<faScalarMatrix> boundedBackwardFaDdtScheme::famDdt
)
);
scalarField limiter(pos(phict) - pos(phict - 1.0));
const scalarField limiter(pos(phict) - pos(phict - 1.0));
scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
scalarField coefft00(limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0)));
scalarField coefft0(coefft + coefft00);
const scalarField coefft(1.0 + limiter*deltaT/(deltaT + deltaT0));
const scalarField coefft00
(
limiter*deltaT*deltaT/(deltaT0*(deltaT + deltaT0))
);
const scalarField coefft0(coefft + coefft00);
fam.diag() = (coefft*rDeltaT)*rho.primitiveField()*mesh().S();

View File

@ -76,10 +76,8 @@ class boundedBackwardFaDdtScheme
{
return GREAT;
}
else
{
return deltaT0_();
}
return deltaT0_();
}
@ -114,7 +112,7 @@ public:
// Member Functions
//- Return mesh reference
const faMesh& mesh() const
const faMesh& mesh() const noexcept
{
return fa::faDdtScheme<scalar>::mesh();
}

View File

@ -28,6 +28,7 @@ License
#include "fa.H"
#include "HashTable.H"
#include "IOobject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -39,6 +40,24 @@ namespace Foam
namespace fa
{
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class Type>
IOobject faDdtScheme<Type>::fieldIOobject(const word& fldName) const
{
const auto& pMesh = mesh().mesh();
return IOobject
(
fldName,
pMesh.time().timeName(),
pMesh,
IOobject::NO_READ,
IOobject::NO_WRITE
);
}
// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
template<class Type>

View File

@ -49,10 +49,12 @@ SourceFiles
namespace Foam
{
// Forward Declarations
template<class Type>
class faMatrix;
class faMesh;
class IOobject;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -68,15 +70,18 @@ class faDdtScheme
:
public refCount
{
protected:
// Protected data
// Protected Data
//- Reference to the mesh
const faMesh& mesh_;
// Private Member Functions
// Protected Member Functions
//- Helper to construct IOobject for field and current time
IOobject fieldIOobject(const word& fldName) const;
//- No copy construct
faDdtScheme(const faDdtScheme&) = delete;
@ -135,10 +140,7 @@ public:
// Member Functions
//- Return mesh reference
const faMesh& mesh() const
{
return mesh_;
}
const faMesh& mesh() const noexcept { return mesh_; }
virtual tmp<GeometricField<Type, faPatchField, areaMesh>> facDdt
(

View File

@ -50,12 +50,7 @@ steadyStateFaDdtScheme<Type>::facDdt
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
IOobject
(
"ddt("+dt.name()+')',
mesh()().time().timeName(),
mesh()()
),
this->fieldIOobject("ddt("+dt.name()+')'),
mesh(),
dimensioned<Type>(dt.dimensions()/dimTime, Zero)
);
@ -71,12 +66,7 @@ steadyStateFaDdtScheme<Type>::facDdt0
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
IOobject
(
"ddt("+dt.name()+')',
mesh()().time().timeName(),
mesh()()
),
this->fieldIOobject("ddt("+dt.name()+')'),
mesh(),
dimensioned<Type>(dt.dimensions()/dimTime, Zero)
);
@ -92,12 +82,7 @@ steadyStateFaDdtScheme<Type>::facDdt
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
IOobject
(
"ddt("+vf.name()+')',
mesh()().time().timeName(),
mesh()()
),
this->fieldIOobject("ddt("+vf.name()+')'),
mesh(),
dimensioned<Type>(vf.dimensions()/dimTime, Zero)
);
@ -113,12 +98,7 @@ steadyStateFaDdtScheme<Type>::facDdt0
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
IOobject
(
"ddt0("+vf.name()+')',
mesh()().time().timeName(),
mesh()()
),
this->fieldIOobject("ddt0("+vf.name()+')'),
mesh(),
dimensioned<Type>(vf.dimensions()/dimTime, Zero)
);
@ -135,17 +115,13 @@ steadyStateFaDdtScheme<Type>::facDdt
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
IOobject
(
"ddt("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()()
),
this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')'),
mesh(),
dimensioned<Type>(rho.dimensions()*vf.dimensions()/dimTime, Zero)
);
}
template<class Type>
tmp<GeometricField<Type, faPatchField, areaMesh>>
steadyStateFaDdtScheme<Type>::facDdt0
@ -156,12 +132,7 @@ steadyStateFaDdtScheme<Type>::facDdt0
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
IOobject
(
"ddt0("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()()
),
this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')'),
mesh(),
dimensioned<Type>(rho.dimensions()*vf.dimensions()/dimTime, Zero)
);
@ -178,12 +149,7 @@ steadyStateFaDdtScheme<Type>::facDdt
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
IOobject
(
"ddt("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()()
),
this->fieldIOobject("ddt("+rho.name()+','+vf.name()+')'),
mesh(),
dimensioned<Type>(rho.dimensions()*vf.dimensions()/dimTime, Zero)
);
@ -200,17 +166,13 @@ steadyStateFaDdtScheme<Type>::facDdt0
{
return tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
IOobject
(
"ddt0("+rho.name()+','+vf.name()+')',
mesh()().time().timeName(),
mesh()()
),
this->fieldIOobject("ddt0("+rho.name()+','+vf.name()+')'),
mesh(),
dimensioned<Type>(rho.dimensions()*vf.dimensions()/dimTime, Zero)
);
}
template<class Type>
tmp<faMatrix<Type>>
steadyStateFaDdtScheme<Type>::famDdt
@ -218,16 +180,11 @@ steadyStateFaDdtScheme<Type>::famDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
return tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
vf.dimensions()*dimArea/dimTime
)
vf,
vf.dimensions()*dimArea/dimTime
);
return tfam;
}
@ -239,16 +196,11 @@ steadyStateFaDdtScheme<Type>::famDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
return tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
)
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
);
return tfam;
}
@ -260,16 +212,11 @@ steadyStateFaDdtScheme<Type>::famDdt
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
tmp<faMatrix<Type>> tfam
return tmp<faMatrix<Type>>::New
(
new faMatrix<Type>
(
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
)
vf,
rho.dimensions()*vf.dimensions()*dimArea/dimTime
);
return tfam;
}

View File

@ -50,6 +50,7 @@ SourceFiles
namespace Foam
{
// Forward Declarations
template<class Type>
class faMatrix;
@ -69,16 +70,18 @@ class divScheme
:
public refCount
{
protected:
// Protected data
//- Reference to the mesh
const faMesh& mesh_;
//- Edge-interpolation scheme
tmp<edgeInterpolationScheme<Type>> tinterpScheme_;
// Private Member Functions
// Protected Member Functions
//- No copy construct
divScheme(const divScheme&) = delete;
@ -135,10 +138,7 @@ public:
// Member Functions
//- Return mesh reference
const faMesh& mesh() const
{
return mesh_;
}
const faMesh& mesh() const noexcept { return mesh_; }
virtual tmp
<

View File

@ -50,23 +50,19 @@ average
{
const faMesh& mesh = ssf.mesh();
tmp<GeometricField<Type, faPatchField, areaMesh>> taverage
auto taverage = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
IOobject
(
IOobject
(
"average("+ssf.name()+')',
ssf.instance(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
ssf.dimensions()
)
"average("+ssf.name()+')',
ssf.instance(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
ssf.dimensions()
);
GeometricField<Type, faPatchField, areaMesh>& av = taverage.ref();
av.ref() =

View File

@ -34,7 +34,6 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef facAverage_H
#define facAverage_H

View File

@ -37,7 +37,6 @@ Author
\*---------------------------------------------------------------------------*/
#ifndef facDdt_H
#define facDdt_H

View File

@ -34,7 +34,6 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef facDiv_H
#define facDiv_H

View File

@ -50,20 +50,17 @@ edgeIntegrate
{
const faMesh& mesh = ssf.mesh();
tmp<GeometricField<Type, faPatchField, areaMesh>> tvf
auto tvf = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
IOobject
(
IOobject
(
"edgeIntegrate("+ssf.name()+')',
ssf.instance(),
ssf.db()
),
mesh,
dimensioned<Type>(ssf.dimensions()/dimArea, Zero),
zeroGradientFaPatchField<Type>::typeName
)
"edgeIntegrate("+ssf.name()+')',
ssf.instance(),
ssf.db()
),
mesh,
dimensioned<Type>(ssf.dimensions()/dimArea, Zero),
zeroGradientFaPatchField<Type>::typeName
);
GeometricField<Type, faPatchField, areaMesh>& vf = tvf.ref();
@ -122,20 +119,17 @@ edgeSum
{
const faMesh& mesh = ssf.mesh();
tmp<GeometricField<Type, faPatchField, areaMesh>> tvf
auto tvf = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
(
new GeometricField<Type, faPatchField, areaMesh>
IOobject
(
IOobject
(
"edgeSum("+ssf.name()+')',
ssf.instance(),
ssf.db()
),
mesh,
dimensioned<Type>(ssf.dimensions(), Zero),
zeroGradientFaPatchField<Type>::typeName
)
"edgeSum("+ssf.name()+')',
ssf.instance(),
ssf.db()
),
mesh,
dimensioned<Type>(ssf.dimensions(), Zero),
zeroGradientFaPatchField<Type>::typeName
);
GeometricField<Type, faPatchField, areaMesh>& vf = tvf.ref();

View File

@ -35,7 +35,6 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef facEdgeIntegrate_H
#define facEdgeIntegrate_H

View File

@ -34,7 +34,6 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef facGrad_H
#define facGrad_H

View File

@ -34,7 +34,6 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef facLaplacian_H
#define facLaplacian_H

View File

@ -34,7 +34,6 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef facLnGrad_H
#define facLnGrad_H

View File

@ -49,6 +49,7 @@ SourceFiles
namespace Foam
{
// Forward Declarations
class faMesh;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -65,8 +66,9 @@ class gradScheme
:
public refCount
{
// Private data
// Private Data
//- Reference to the mesh
const faMesh& mesh_;
@ -119,10 +121,7 @@ public:
// Member Functions
//- Return mesh reference
const faMesh& mesh() const
{
return mesh_;
}
const faMesh& mesh() const noexcept { return mesh_; }
//- Calculate and return the grad of the given field
virtual tmp

View File

@ -62,8 +62,9 @@ class gaussGrad
:
public fa::gradScheme<Type>
{
// Private data
// Private Data
//- Edge-interpolation scheme
tmp<edgeInterpolationScheme<Type>> tinterpScheme_;

View File

@ -63,22 +63,19 @@ leastSquaresFaGrad<Type>::grad
const faMesh& mesh = vsf.mesh();
tmp<GeometricField<GradType, faPatchField, areaMesh>> tlsGrad
auto tlsGrad = tmp<GeometricField<GradType, faPatchField, areaMesh>>::New
(
new GeometricField<GradType, faPatchField, areaMesh>
IOobject
(
IOobject
(
"grad(" + vsf.name() + ')',
vsf.instance(),
vsf.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensioned<GradType>(vsf.dimensions()/dimLength, Zero),
zeroGradientFaPatchField<GradType>::typeName
)
"grad(" + vsf.name() + ')',
vsf.instance(),
vsf.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensioned<GradType>(vsf.dimensions()/dimLength, Zero),
zeroGradientFaPatchField<GradType>::typeName
);
GeometricField<GradType, faPatchField, areaMesh>& lsGrad = tlsGrad.ref();
@ -93,10 +90,10 @@ leastSquaresFaGrad<Type>::grad
forAll(own, edgei)
{
label ownEdgeI = own[edgei];
label neiEdgeI = nei[edgei];
const label ownEdgeI = own[edgei];
const label neiEdgeI = nei[edgei];
Type deltaVsf = vsf[neiEdgeI] - vsf[ownEdgeI];
const Type deltaVsf(vsf[neiEdgeI] - vsf[ownEdgeI]);
lsGrad[ownEdgeI] += ownLs[edgei]*deltaVsf;
lsGrad[neiEdgeI] -= neiLs[edgei]*deltaVsf;
@ -112,7 +109,7 @@ leastSquaresFaGrad<Type>::grad
if (vsf.boundaryField()[patchi].coupled())
{
Field<Type> neiVsf
const Field<Type> neiVsf
(
vsf.boundaryField()[patchi].patchNeighbourField()
);

View File

@ -54,8 +54,8 @@ Foam::leastSquaresFaVectors::leastSquaresFaVectors(const faMesh& mesh)
Foam::leastSquaresFaVectors::~leastSquaresFaVectors()
{
deleteDemandDrivenData(pVectorsPtr_);
deleteDemandDrivenData(nVectorsPtr_);
pVectorsPtr_.reset(nullptr);
nVectorsPtr_.reset(nullptr);
}
@ -66,35 +66,41 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
DebugInFunction
<< "Constructing finite area least square gradient vectors" << nl;
pVectorsPtr_ = new edgeVectorField
pVectorsPtr_.reset
(
IOobject
new edgeVectorField
(
"LeastSquaresP",
mesh().pointsInstance(),
mesh().thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh(),
dimensionedVector(dimless/dimLength, Zero)
IOobject
(
"LeastSquaresP",
mesh().pointsInstance(),
mesh().thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh(),
dimensionedVector(dimless/dimLength, Zero)
)
);
edgeVectorField& lsP = *pVectorsPtr_;
nVectorsPtr_ = new edgeVectorField
nVectorsPtr_.reset
(
IOobject
new edgeVectorField
(
"LeastSquaresN",
mesh().pointsInstance(),
mesh().thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh(),
dimensionedVector(dimless/dimLength, Zero)
IOobject
(
"LeastSquaresN",
mesh().pointsInstance(),
mesh().thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh(),
dimensionedVector(dimless/dimLength, Zero)
)
);
edgeVectorField& lsN = *nVectorsPtr_;
@ -111,11 +117,18 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
const label own = owner[facei];
const label nei = neighbour[facei];
vector d = C[nei] - C[own];
symmTensor wdd = (1.0/magSqr(d))*sqr(d);
vector d(C[nei] - C[own]);
// Do not allow any mag(val) < SMALL
if (mag(d) < SMALL)
{
d = vector::uniform(SMALL);
}
const symmTensor wdd((1.0/magSqr(d))*sqr(d));
dd[own] += wdd;
dd[nei] += wdd;
@ -135,41 +148,34 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
// cause problems with fixedGradient. HJ, 4/Oct/2010
const vectorField pd(p.delta());
if (p.coupled())
forAll(pd, patchFacei)
{
forAll(pd, patchFacei)
{
const vector& d = pd[patchFacei];
const vector& d = pd[patchFacei];
dd[edgeFaces[patchFacei]] +=
(1.0/magSqr(d))*sqr(d);
}
}
else
{
forAll(pd, patchFacei)
{
const vector& d = pd[patchFacei];
dd[edgeFaces[patchFacei]] +=
(1.0/magSqr(d))*sqr(d);
}
dd[edgeFaces[patchFacei]] += (1.0/magSqr(d))*sqr(d);
}
}
// Invert the dd tensor
const symmTensorField invDd(pinv(dd));
const symmTensorField invDd(inv(dd));
// Revisit all faces and calculate the lsP and lsN vectors
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
const label own = owner[facei];
const label nei = neighbour[facei];
vector d = C[nei] - C[own];
scalar magSfByMagSqrd = 1.0/magSqr(d);
vector d(C[nei] - C[own]);
// Do not allow any mag(val) < SMALL
if (mag(d) < SMALL)
{
d = vector::uniform(SMALL);
}
const scalar magSfByMagSqrd = 1.0/magSqr(d);
lsP[facei] = magSfByMagSqrd*(invDd[own] & d);
lsN[facei] = -magSfByMagSqrd*(invDd[nei] & d);
@ -187,27 +193,12 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
// Build the d-vectors
const vectorField pd(p.delta());
if (p.coupled())
forAll(pd, patchFacei)
{
forAll(pd, patchFacei)
{
const vector& d = pd[patchFacei];
const vector& d = pd[patchFacei];
patchLsP[patchFacei] =
(1.0/magSqr(d))
*(invDd[edgeFaces[patchFacei]] & d);
}
}
else
{
forAll(pd, patchFacei)
{
const vector& d = pd[patchFacei];
patchLsP[patchFacei] =
(1.0/magSqr(d))
*(invDd[edgeFaces[patchFacei]] & d);
}
patchLsP[patchFacei] =
(1.0/magSqr(d))*(invDd[edgeFaces[patchFacei]] & d);
}
}
@ -243,8 +234,8 @@ bool Foam::leastSquaresFaVectors::movePoints()
DebugInFunction
<< "Clearing least square data" << nl;
deleteDemandDrivenData(pVectorsPtr_);
deleteDemandDrivenData(nVectorsPtr_);
pVectorsPtr_.reset(nullptr);
nVectorsPtr_.reset(nullptr);
return true;
}

View File

@ -61,8 +61,8 @@ class leastSquaresFaVectors
// Private data
//- Least-squares gradient vectors
mutable edgeVectorField* pVectorsPtr_;
mutable edgeVectorField* nVectorsPtr_;
mutable std::unique_ptr<edgeVectorField> pVectorsPtr_;
mutable std::unique_ptr<edgeVectorField> nVectorsPtr_;
// Private member functions

View File

@ -65,6 +65,7 @@ class edgeLimitedGrad
{
// Private Data
//- Basic gradient scheme
tmp<fa::gradScheme<Type>> basicGradScheme_;
//- Limiter coefficient

View File

@ -97,19 +97,19 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::grad
// create limiter
scalarField limiter(vsf.internalField().size(), 1.0);
scalar rk = (1.0/k_ - 1.0);
const scalar rk = (1.0/k_ - 1.0);
forAll(owner, edgei)
{
label own = owner[edgei];
label nei = neighbour[edgei];
const label own = owner[edgei];
const label nei = neighbour[edgei];
scalar vsfOwn = vsf[own];
scalar vsfNei = vsf[nei];
const scalar vsfOwn = vsf[own];
const scalar vsfNei = vsf[nei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
@ -145,14 +145,14 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::grad
forAll(pOwner, pEdgei)
{
label own = pOwner[pEdgei];
const label own = pOwner[pEdgei];
scalar vsfOwn = vsf[own];
scalar vsfNei = psfNei[pEdgei];
const scalar vsfOwn = vsf[own];
const scalar vsfNei = psfNei[pEdgei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
@ -168,14 +168,14 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::grad
{
forAll(pOwner, pEdgei)
{
label own = pOwner[pEdgei];
const label own = pOwner[pEdgei];
scalar vsfOwn = vsf[own];
scalar vsfNei = psf[pEdgei];
const scalar vsfOwn = vsf[own];
const scalar vsfNei = psf[pEdgei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
@ -231,25 +231,25 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::grad
// create limiter
scalarField limiter(vvf.internalField().size(), 1.0);
scalar rk = (1.0/k_ - 1.0);
const scalar rk = (1.0/k_ - 1.0);
forAll(owner, edgei)
{
label own = owner[edgei];
label nei = neighbour[edgei];
const label own = owner[edgei];
const label nei = neighbour[edgei];
vector vvfOwn = vvf[own];
vector vvfNei = vvf[nei];
const vector vvfOwn(vvf[own]);
const vector vvfNei(vvf[nei]);
// owner side
vector gradf = (Cf[edgei] - C[own]) & g[own];
vector gradf((Cf[edgei] - C[own]) & g[own]);
scalar vsfOwn = gradf & vvfOwn;
scalar vsfNei = gradf & vvfNei;
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
@ -294,19 +294,19 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::grad
forAll(pOwner, pEdgei)
{
label own = pOwner[pEdgei];
const label own = pOwner[pEdgei];
vector vvfOwn = vvf[own];
vector vvfNei = psfNei[pEdgei];
const vector vvfOwn(vvf[own]);
const vector vvfNei(psfNei[pEdgei]);
vector gradf = (pCf[pEdgei] - C[own]) & g[own];
const vector gradf((pCf[pEdgei] - C[own]) & g[own]);
scalar vsfOwn = gradf & vvfOwn;
scalar vsfNei = gradf & vvfNei;
const scalar vsfOwn = gradf & vvfOwn;
const scalar vsfNei = gradf & vvfNei;
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
@ -322,19 +322,19 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::grad
{
forAll(pOwner, pEdgei)
{
label own = pOwner[pEdgei];
const label own = pOwner[pEdgei];
vector vvfOwn = vvf[own];
vector vvfNei = psf[pEdgei];
const vector vvfOwn(vvf[own]);
const vector vvfNei(psf[pEdgei]);
vector gradf = (pCf[pEdgei] - C[own]) & g[own];
const vector gradf((pCf[pEdgei] - C[own]) & g[own]);
scalar vsfOwn = gradf & vvfOwn;
scalar vsfNei = gradf & vvfNei;
const scalar vsfOwn = gradf & vvfOwn;
const scalar vsfNei = gradf & vvfNei;
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;

View File

@ -68,6 +68,7 @@ class faceLimitedGrad
{
// Private Data
//- Basic gradient scheme
tmp<fa::gradScheme<Type>> basicGradScheme_;
//- Limiter coefficient

View File

@ -121,11 +121,11 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
const label own = owner[facei];
const label nei = neighbour[facei];
scalar vsfOwn = vsf[own];
scalar vsfNei = vsf[nei];
const scalar vsfOwn = vsf[own];
const scalar vsfNei = vsf[nei];
maxVsf[own] = max(maxVsf[own], vsfNei);
minVsf[own] = min(minVsf[own], vsfNei);
@ -149,8 +149,8 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
scalar vsfNei = psfNei[pFacei];
const label own = pOwner[pFacei];
const scalar vsfNei = psfNei[pFacei];
maxVsf[own] = max(maxVsf[own], vsfNei);
minVsf[own] = min(minVsf[own], vsfNei);
@ -160,8 +160,8 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
{
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
scalar vsfNei = psf[pFacei];
const label own = pOwner[pFacei];
const scalar vsfNei = psf[pFacei];
maxVsf[own] = max(maxVsf[own], vsfNei);
minVsf[own] = min(minVsf[own], vsfNei);
@ -174,7 +174,7 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
if (k_ < 1.0)
{
scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
const scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
maxVsf += maxMinVsf;
minVsf -= maxMinVsf;
}
@ -185,8 +185,8 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
const label own = owner[facei];
const label nei = neighbour[facei];
// owner side
limitEdge
@ -214,7 +214,7 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::grad
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
const label own = pOwner[pFacei];
limitEdge
(
@ -270,8 +270,8 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
const label own = owner[facei];
const label nei = neighbour[facei];
const vector& vsfOwn = vsf[own];
const vector& vsfNei = vsf[nei];
@ -293,11 +293,11 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
if (psf.coupled())
{
vectorField psfNei(psf.patchNeighbourField());
const vectorField psfNei(psf.patchNeighbourField());
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
const label own = pOwner[pFacei];
const vector& vsfNei = psfNei[pFacei];
maxVsf[own] = max(maxVsf[own], vsfNei);
@ -308,7 +308,7 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
{
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
const label own = pOwner[pFacei];
const vector& vsfNei = psf[pFacei];
maxVsf[own] = max(maxVsf[own], vsfNei);
@ -322,7 +322,7 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
if (k_ < 1.0)
{
vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
const vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
maxVsf += maxMinVsf;
minVsf -= maxMinVsf;
@ -336,8 +336,8 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
const label own = owner[facei];
const label nei = neighbour[facei];
// owner side
limitEdge
@ -365,7 +365,7 @@ tmp<areaTensorField> faceLimitedGrad<vector>::grad
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
const label own = pOwner[pFacei];
limitEdge
(

View File

@ -59,19 +59,16 @@ correctedLnGrad<Type>::correction
{
const faMesh& mesh = this->mesh();
tmp<GeometricField<Type, faePatchField, edgeMesh>> tssf
auto tssf = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
(
new GeometricField<Type, faePatchField, edgeMesh>
IOobject
(
IOobject
(
"lnGradCorr("+vf.name()+')',
vf.instance(),
vf.db()
),
mesh,
vf.dimensions()*mesh.deltaCoeffs().dimensions()
)
"lnGradCorr("+vf.name()+')',
vf.instance(),
vf.db()
),
mesh,
vf.dimensions()*mesh.deltaCoeffs().dimensions()
);
GeometricField<Type, faePatchField, edgeMesh>& ssf = tssf.ref();
@ -80,7 +77,7 @@ correctedLnGrad<Type>::correction
ssf.replace
(
cmpt,
mesh.correctionVectors()
mesh.nonOrthCorrectionVector()
& linearEdgeInterpolation
<
typename

View File

@ -102,7 +102,7 @@ public:
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
virtual bool corrected() const noexcept
{
return !this->mesh().orthogonal();
}

View File

@ -60,23 +60,20 @@ fourthLnGrad<Type>::correction
{
const faMesh& mesh = this->mesh();
tmp<GeometricField<Type, faePatchField, edgeMesh>> tcorr
auto tcorr = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
(
new GeometricField<Type, faePatchField, edgeMesh>
IOobject
(
IOobject
(
"lnGradCorr("+vf.name()+')',
vf.instance(),
vf.db()
),
mesh,
vf.dimensions()*this->mesh().deltaCoeffs().dimensions()
)
"lnGradCorr("+vf.name()+')',
vf.instance(),
vf.db()
),
mesh,
vf.dimensions()*this->mesh().deltaCoeffs().dimensions()
);
GeometricField<Type, faePatchField, edgeMesh>& corr = tcorr.ref();
edgeVectorField m(mesh.Le()/mesh.magLe());
const edgeVectorField m(mesh.Le()/mesh.magLe());
for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; ++cmpt)
{

View File

@ -102,10 +102,7 @@ public:
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return true;
}
virtual bool corrected() const noexcept { return true; }
//- Return the explicit correction to the fourthLnGrad
// for the given field

View File

@ -70,17 +70,17 @@ limitedLnGrad<Type>::correction
*mag(lnGradScheme<Type>::lnGrad(vf, deltaCoeffs(vf), "orthSnGrad"))
/(
(1 - limitCoeff_)*mag(corr)
+ dimensionedScalar("small", corr.dimensions(), SMALL)
+ dimensionedScalar(corr.dimensions(), SMALL)
),
dimensionedScalar("one", dimless, 1.0)
dimensionedScalar(dimless, 1.0)
)
);
if (fa::debug)
{
Info<< "limitedLnGrad :: limiter min: " << min(limiter.internalField())
<< " max: "<< max(limiter.internalField())
<< " avg: " << average(limiter.internalField()) << endl;
Info<< "limitedLnGrad :: limiter min: " << gMin(limiter.internalField())
<< " max: "<< gMax(limiter.internalField())
<< " avg: " << gAverage(limiter.internalField()) << endl;
}
return limiter*corr;

View File

@ -68,7 +68,7 @@ class limitedLnGrad
:
public lnGradScheme<Type>
{
// Private data
// Private Data
//- Limiter. 0 = no limiting, 1 = full limiting
scalar limitCoeff_;
@ -128,7 +128,7 @@ public:
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
virtual bool corrected() const noexcept
{
return !this->mesh().orthogonal();
}

View File

@ -106,22 +106,18 @@ lnGradScheme<Type>::lnGrad
{
const faMesh& mesh = vf.mesh();
// construct GeometricField<Type, faePatchField, edgeMesh>
tmp<GeometricField<Type, faePatchField, edgeMesh>> tssf
auto tssf = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
(
new GeometricField<Type, faePatchField, edgeMesh>
IOobject
(
IOobject
(
lnGradName + "("+vf.name()+')',
vf.instance(),
vf.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
vf.dimensions()*tdeltaCoeffs().dimensions()
)
lnGradName + "("+vf.name()+')',
vf.instance(),
vf.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
vf.dimensions()*tdeltaCoeffs().dimensions()
);
GeometricField<Type, faePatchField, edgeMesh>& ssf = tssf.ref();
@ -132,17 +128,17 @@ lnGradScheme<Type>::lnGrad
const labelUList& owner = mesh.owner();
const labelUList& neighbour = mesh.neighbour();
forAll(owner, faceI)
forAll(owner, facei)
{
ssf[faceI] =
deltaCoeffs[faceI]*(vf[neighbour[faceI]] - vf[owner[faceI]]);
ssf[facei] =
deltaCoeffs[facei]*(vf[neighbour[facei]] - vf[owner[facei]]);
}
auto& ssfb = ssf.boundaryFieldRef();
forAll(vf.boundaryField(), patchI)
forAll(vf.boundaryField(), patchi)
{
ssfb[patchI] = vf.boundaryField()[patchI].snGrad();
ssfb[patchi] = vf.boundaryField()[patchi].snGrad();
}
return tssf;

View File

@ -48,6 +48,7 @@ SourceFiles
namespace Foam
{
// Forward Declarations
class faMesh;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -64,7 +65,7 @@ class lnGradScheme
:
public refCount
{
// Private data
// Private Data
//- Hold reference to mesh
const faMesh& mesh_;
@ -116,10 +117,7 @@ public:
// Member Functions
//- Return mesh reference
const faMesh& mesh() const
{
return mesh_;
}
const faMesh& mesh() const noexcept { return mesh_; }
//- Return the lnGrad of the given cell field
// with the given weighting factors
@ -138,10 +136,7 @@ public:
) const = 0;
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return false;
}
virtual bool corrected() const noexcept { return false; }
//- Return the explicit correction to the lnGrad
// for the given field

View File

@ -101,10 +101,7 @@ public:
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return false;
}
virtual bool corrected() const noexcept { return false; }
//- Return the explicit correction to the uncorrectedLnGrad
//- for the given field

View File

@ -100,14 +100,14 @@ Foam::fac::interpolate
Istream& schemeData
)
{
# ifdef DEBUGInterpolations
#ifdef DEBUGInterpolations
if (edgeInterpolation::debug)
{
InfoInFunction
<< "interpolating GeometricField<Type, faPatchField, areaMesh> "
<< endl;
}
# endif
#endif
return scheme<Type>(faceFlux, schemeData)().interpolate(vf);
}
@ -122,7 +122,7 @@ Foam::fac::interpolate
const word& name
)
{
# ifdef DEBUGInterpolations
#ifdef DEBUGInterpolations
if (edgeInterpolation::debug)
{
InfoInFunction
@ -130,7 +130,7 @@ Foam::fac::interpolate
<< "using " << name
<< endl;
}
# endif
#endif
return scheme<Type>(faceFlux, name)().interpolate(vf);
}
@ -199,14 +199,14 @@ Foam::fac::interpolate
Istream& schemeData
)
{
# ifdef DEBUGInterpolations
#ifdef DEBUGInterpolations
if (edgeInterpolation::debug)
{
InfoInFunction
<< "interpolating GeometricField<Type, faPatchField, areaMesh> "
<< endl;
}
# endif
#endif
return scheme<Type>(vf.mesh(), schemeData)().interpolate(vf);
}
@ -220,7 +220,7 @@ Foam::fac::interpolate
const word& name
)
{
# ifdef DEBUGInterpolations
#ifdef DEBUGInterpolations
if (edgeInterpolation::debug)
{
InfoInFunction
@ -228,7 +228,7 @@ Foam::fac::interpolate
<< "using " << name
<< endl;
}
# endif
#endif
return scheme<Type>(vf.mesh(), name)().interpolate(vf);
}
@ -257,7 +257,7 @@ Foam::fac::interpolate
const GeometricField<Type, faPatchField, areaMesh>& vf
)
{
# ifdef DEBUGInterpolations
#ifdef DEBUGInterpolations
if (edgeInterpolation::debug)
{
InfoInFunction
@ -265,7 +265,7 @@ Foam::fac::interpolate
<< "using run-time selected scheme"
<< endl;
}
# endif
#endif
return interpolate(vf, "interpolate(" + vf.name() + ')');
}

View File

@ -31,6 +31,7 @@ License
#include "edgeFields.H"
#include "demandDrivenData.H"
#include "faPatchFields.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -41,13 +42,13 @@ namespace Foam
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void Foam::edgeInterpolation::clearOut()
void Foam::edgeInterpolation::clearOut() const
{
deleteDemandDrivenData(lPN_);
deleteDemandDrivenData(weightingFactors_);
deleteDemandDrivenData(differenceFactors_);
deleteDemandDrivenData(correctionVectors_);
deleteDemandDrivenData(skewCorrectionVectors_);
lPNptr_.reset(nullptr);
weightsPtr_.reset(nullptr);
deltaCoeffsPtr_.reset(nullptr);
nonOrthCorrectionVectorsPtr_.reset(nullptr);
skewCorrectionVectorsPtr_.reset(nullptr);
}
@ -56,11 +57,11 @@ void Foam::edgeInterpolation::clearOut()
Foam::edgeInterpolation::edgeInterpolation(const faMesh& fam)
:
faMesh_(fam),
lPN_(nullptr),
weightingFactors_(nullptr),
differenceFactors_(nullptr),
correctionVectors_(nullptr),
skewCorrectionVectors_(nullptr),
lPNptr_(nullptr),
weightsPtr_(nullptr),
deltaCoeffsPtr_(nullptr),
nonOrthCorrectionVectorsPtr_(nullptr),
skewCorrectionVectorsPtr_(nullptr),
orthogonal_(false),
skew_(true)
{}
@ -78,40 +79,40 @@ Foam::edgeInterpolation::~edgeInterpolation()
const Foam::edgeScalarField& Foam::edgeInterpolation::lPN() const
{
if (!lPN_)
if (!lPNptr_)
{
makeLPN();
}
return (*lPN_);
return (*lPNptr_);
}
const Foam::edgeScalarField& Foam::edgeInterpolation::weights() const
{
if (!weightingFactors_)
if (!weightsPtr_)
{
makeWeights();
}
return (*weightingFactors_);
return (*weightsPtr_);
}
const Foam::edgeScalarField& Foam::edgeInterpolation::deltaCoeffs() const
{
if (!differenceFactors_)
if (!deltaCoeffsPtr_)
{
makeDeltaCoeffs();
}
return (*differenceFactors_);
return (*deltaCoeffsPtr_);
}
bool Foam::edgeInterpolation::orthogonal() const
{
if (orthogonal_ == false && !correctionVectors_)
if (orthogonal_ == false && !nonOrthCorrectionVectorsPtr_)
{
makeCorrectionVectors();
}
@ -120,7 +121,8 @@ bool Foam::edgeInterpolation::orthogonal() const
}
const Foam::edgeVectorField& Foam::edgeInterpolation::correctionVectors() const
const Foam::edgeVectorField&
Foam::edgeInterpolation::nonOrthCorrectionVector() const
{
if (orthogonal())
{
@ -129,13 +131,13 @@ const Foam::edgeVectorField& Foam::edgeInterpolation::correctionVectors() const
<< abort(FatalError);
}
return (*correctionVectors_);
return (*nonOrthCorrectionVectorsPtr_);
}
bool Foam::edgeInterpolation::skew() const
{
if (skew_ == true && !skewCorrectionVectors_)
if (skew_ == true && !skewCorrectionVectorsPtr_)
{
makeSkewCorrectionVectors();
}
@ -154,48 +156,66 @@ Foam::edgeInterpolation::skewCorrectionVectors() const
<< abort(FatalError);
}
return (*skewCorrectionVectors_);
return (*skewCorrectionVectorsPtr_);
}
bool Foam::edgeInterpolation::movePoints() const
{
deleteDemandDrivenData(lPN_);
deleteDemandDrivenData(weightingFactors_);
deleteDemandDrivenData(differenceFactors_);
orthogonal_ = false;
deleteDemandDrivenData(correctionVectors_);
skew_ = true;
deleteDemandDrivenData(skewCorrectionVectors_);
clearOut();
return true;
}
const Foam::vector& Foam::edgeInterpolation::skewCorr(const label edgei) const
{
#ifdef FA_SKEW_CORRECTION
return
(
skewCorrectionVectorsPtr_
? (*skewCorrectionVectorsPtr_)[edgei]
: pTraits<vector>::zero
);
#else
return (*skewCorrectionVectorsPtr_)[edgei];
#endif
}
void Foam::edgeInterpolation::makeLPN() const
{
DebugInFunction
<< "Constructing geodesic distance between points P and N"
<< endl;
const polyMesh& pMesh = mesh().mesh();
lPN_ = new edgeScalarField
lPNptr_.reset
(
IOobject
new edgeScalarField
(
"lPN",
faMesh_.time().constant(),
faMesh_(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh(),
dimLength
IOobject
(
"lPN",
pMesh.pointsInstance(),
pMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh(),
dimensionedScalar(dimLength, Zero)
)
);
edgeScalarField& lPN = *lPN_;
edgeScalarField& lPN = *lPNptr_;
// Set local references to mesh data
const edgeVectorField& edgeCentres = mesh().edgeCentres();
@ -205,43 +225,47 @@ void Foam::edgeInterpolation::makeLPN() const
scalarField& lPNIn = lPN.primitiveFieldRef();
forAll(owner, edgeI)
{
vector curSkewCorrVec(Zero);
// Calculate skewness correction vectors if necessary
(void) skew();
if (skew())
{
curSkewCorrVec = skewCorrectionVectors()[edgeI];
}
forAll(owner, edgei)
{
const vector& skewCorrEdge = skewCorr(edgei);
scalar lPE =
mag
(
edgeCentres[edgeI]
- curSkewCorrVec
- faceCentres[owner[edgeI]]
edgeCentres[edgei]
- skewCorrEdge
- faceCentres[owner[edgei]]
);
scalar lEN =
mag
(
faceCentres[neighbour[edgeI]]
- edgeCentres[edgeI]
+ curSkewCorrVec
faceCentres[neighbour[edgei]]
- edgeCentres[edgei]
+ skewCorrEdge
);
lPNIn[edgeI] = (lPE + lEN);
lPNIn[edgei] = (lPE + lEN);
// Do not allow any mag(val) < SMALL
if (mag(lPNIn[edgei]) < SMALL)
{
lPNIn[edgei] = SMALL;
}
}
forAll(lPN.boundaryField(), patchI)
forAll(lPN.boundaryField(), patchi)
{
mesh().boundary()[patchI].makeDeltaCoeffs
mesh().boundary()[patchi].makeDeltaCoeffs
(
lPN.boundaryFieldRef()[patchI]
lPN.boundaryFieldRef()[patchi]
);
lPN.boundaryFieldRef()[patchI] = 1.0/lPN.boundaryField()[patchI];
lPN.boundaryFieldRef()[patchi] = 1.0/lPN.boundaryField()[patchi];
}
@ -257,22 +281,26 @@ void Foam::edgeInterpolation::makeWeights() const
<< "Constructing weighting factors for edge interpolation"
<< endl;
const polyMesh& pMesh = mesh().mesh();
weightingFactors_ = new edgeScalarField
weightsPtr_.reset
(
IOobject
new edgeScalarField
(
"weightingFactors",
mesh()().pointsInstance(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh(),
dimless
IOobject
(
"weights",
pMesh.pointsInstance(),
pMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh(),
dimensionedScalar(dimless, 1)
)
);
edgeScalarField& weightingFactors = *weightingFactors_;
edgeScalarField& weights = *weightsPtr_;
// Set local references to mesh data
@ -281,49 +309,44 @@ void Foam::edgeInterpolation::makeWeights() const
const labelUList& owner = mesh().owner();
const labelUList& neighbour = mesh().neighbour();
scalarField& weightingFactorsIn = weightingFactors.primitiveFieldRef();
scalarField& weightsIn = weights.primitiveFieldRef();
forAll(owner, edgeI)
// Calculate skewness correction vectors if necessary
(void) skew();
forAll(owner, edgei)
{
vector curSkewCorrVec(Zero);
if (skew())
{
curSkewCorrVec = skewCorrectionVectors()[edgeI];
}
const vector& skewCorrEdge = skewCorr(edgei);
scalar lPE =
mag
(
edgeCentres[edgeI]
- curSkewCorrVec
- faceCentres[owner[edgeI]]
edgeCentres[edgei]
- skewCorrEdge
- faceCentres[owner[edgei]]
);
scalar lEN =
mag
(
faceCentres[neighbour[edgeI]]
- edgeCentres[edgeI]
+ curSkewCorrVec
faceCentres[neighbour[edgei]]
- edgeCentres[edgei]
+ skewCorrEdge
);
weightingFactorsIn[edgeI] =
lEN
/(
lPE
# ifdef BAD_MESH_STABILISATION
+ VSMALL
# endif
+ lEN
);
// weight = (0,1]
const scalar lPN = lPE + lEN;
if (mag(lPN) > SMALL)
{
weightsIn[edgei] = lEN/lPN;
}
}
forAll(mesh().boundary(), patchI)
forAll(mesh().boundary(), patchi)
{
mesh().boundary()[patchI].makeWeights
mesh().boundary()[patchi].makeWeights
(
weightingFactors.boundaryFieldRef()[patchI]
weights.boundaryFieldRef()[patchi]
);
}
@ -343,22 +366,27 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
// needed to make sure deltaCoeffs are calculated for parallel runs.
weights();
differenceFactors_ = new edgeScalarField
const polyMesh& pMesh = mesh().mesh();
deltaCoeffsPtr_.reset
(
IOobject
new edgeScalarField
(
"differenceFactors_",
mesh()().pointsInstance(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh(),
dimless/dimLength
IOobject
(
"deltaCoeffs",
pMesh.pointsInstance(),
pMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh(),
dimensionedScalar(dimless/dimLength, SMALL)
)
);
edgeScalarField& DeltaCoeffs = *differenceFactors_;
scalarField& dc = DeltaCoeffs.primitiveFieldRef();
edgeScalarField& deltaCoeffs = *deltaCoeffsPtr_;
scalarField& dc = deltaCoeffs.primitiveFieldRef();
// Set local references to mesh data
const edgeVectorField& edgeCentres = mesh().edgeCentres();
@ -370,65 +398,62 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
const edgeList& edges = mesh().edges();
const pointField& points = mesh().points();
// Calculate skewness correction vectors if necessary
(void) skew();
forAll(owner, edgeI)
forAll(owner, edgei)
{
// Edge normal - area normal
vector edgeNormal =
normalised(lengths[edgeI] ^ edges[edgeI].vec(points));
normalised(lengths[edgei] ^ edges[edgei].vec(points));
// Unit delta vector
vector unitDelta =
faceCentres[neighbour[edgeI]]
- faceCentres[owner[edgeI]];
faceCentres[neighbour[edgei]]
- faceCentres[owner[edgei]];
unitDelta.removeCollinear(edgeNormal);
unitDelta.normalise();
// Calc PN arc length
vector curSkewCorrVec(Zero);
if (skew())
{
curSkewCorrVec = skewCorrectionVectors()[edgeI];
}
const vector& skewCorrEdge = skewCorr(edgei);
scalar lPE =
mag
(
edgeCentres[edgeI]
- curSkewCorrVec
- faceCentres[owner[edgeI]]
edgeCentres[edgei]
- skewCorrEdge
- faceCentres[owner[edgei]]
);
scalar lEN =
mag
(
faceCentres[neighbour[edgeI]]
- edgeCentres[edgeI]
+ curSkewCorrVec
faceCentres[neighbour[edgei]]
- edgeCentres[edgei]
+ skewCorrEdge
);
scalar lPN = lPE + lEN;
// Edge normal - area tangent
edgeNormal = normalised(lengths[edgeI]);
edgeNormal = normalised(lengths[edgei]);
// Delta coefficient as per definition
// dc[edgeI] = 1.0/(lPN*(unitDelta & edgeNormal));
// Stabilised form for bad meshes. HJ, 23/Jul/2009
dc[edgeI] = 1.0/max((lPN*(unitDelta & edgeNormal)), 0.05*lPN);
// Do not allow any mag(val) < SMALL
const scalar alpha = lPN*(unitDelta & edgeNormal);
if (mag(alpha) > SMALL)
{
dc[edgei] = scalar(1)/max(alpha, 0.05*lPN);
}
}
forAll(DeltaCoeffs.boundaryField(), patchI)
forAll(deltaCoeffs.boundaryField(), patchi)
{
mesh().boundary()[patchI].makeDeltaCoeffs
mesh().boundary()[patchi].makeDeltaCoeffs
(
DeltaCoeffs.boundaryFieldRef()[patchI]
deltaCoeffs.boundaryFieldRef()[patchi]
);
}
}
@ -440,21 +465,26 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
<< "Constructing non-orthogonal correction vectors"
<< endl;
correctionVectors_ = new edgeVectorField
const polyMesh& pMesh = mesh().mesh();
nonOrthCorrectionVectorsPtr_.reset
(
IOobject
new edgeVectorField
(
"correctionVectors",
mesh()().pointsInstance(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh(),
dimless
IOobject
(
"nonOrthCorrectionVectors",
pMesh.pointsInstance(),
pMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh(),
dimensionedVector(dimless, Zero)
)
);
edgeVectorField& CorrVecs = *correctionVectors_;
edgeVectorField& CorrVecs = *nonOrthCorrectionVectorsPtr_;
// Set local references to mesh data
const areaVectorField& faceCentres = mesh().areaCentres();
@ -467,74 +497,78 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
const edgeList& edges = mesh().edges();
const pointField& points = mesh().points();
scalarField deltaCoeffs(owner.size());
scalarField deltaCoeffs(owner.size(), SMALL);
vectorField& CorrVecsIn = CorrVecs.primitiveFieldRef();
forAll(owner, edgeI)
forAll(owner, edgei)
{
// Edge normal - area normal
vector edgeNormal =
normalised(lengths[edgeI] ^ edges[edgeI].vec(points));
normalised(lengths[edgei] ^ edges[edgei].vec(points));
// Unit delta vector
vector unitDelta =
faceCentres[neighbour[edgeI]]
- faceCentres[owner[edgeI]];
faceCentres[neighbour[edgei]]
- faceCentres[owner[edgei]];
unitDelta.removeCollinear(edgeNormal);
unitDelta.normalise();
// Edge normal - area tangent
edgeNormal = normalised(lengths[edgeI]);
edgeNormal = normalised(lengths[edgei]);
// Delta coeffs
deltaCoeffs[edgeI] = 1.0/(unitDelta & edgeNormal);
// Do not allow any mag(val) < SMALL
const scalar alpha = unitDelta & edgeNormal;
if (mag(alpha) > SMALL)
{
deltaCoeffs[edgei] = scalar(1)/alpha;
}
// Edge correction vector
CorrVecsIn[edgeI] =
CorrVecsIn[edgei] =
edgeNormal
- deltaCoeffs[edgeI]*unitDelta;
- deltaCoeffs[edgei]*unitDelta;
}
edgeVectorField::Boundary& CorrVecsbf = CorrVecs.boundaryFieldRef();
forAll(CorrVecs.boundaryField(), patchI)
forAll(CorrVecs.boundaryField(), patchi)
{
mesh().boundary()[patchI].makeCorrectionVectors(CorrVecsbf[patchI]);
mesh().boundary()[patchi].makeCorrectionVectors(CorrVecsbf[patchi]);
}
scalar NonOrthogCoeff = 0.0;
constexpr scalar maxNonOrthoRatio = 0.1;
scalar nonOrthoCoeff = 0;
if (owner.size() > 0)
{
scalarField sinAlpha(deltaCoeffs*mag(CorrVecs.internalField()));
forAll(sinAlpha, edgeI)
forAll(sinAlpha, edgei)
{
sinAlpha[edgeI] = max(-1, min(sinAlpha[edgeI], 1));
sinAlpha[edgei] = max(-1, min(sinAlpha[edgei], 1));
}
NonOrthogCoeff = max(Foam::asin(sinAlpha)*180.0/M_PI);
nonOrthoCoeff = max(Foam::asin(sinAlpha)*radToDeg());
}
reduce(NonOrthogCoeff, maxOp<scalar>());
reduce(nonOrthoCoeff, maxOp<scalar>());
DebugInFunction
<< "non-orthogonality coefficient = " << NonOrthogCoeff << " deg."
<< "non-orthogonality coefficient = " << nonOrthoCoeff << " deg."
<< endl;
if (NonOrthogCoeff < 0.1)
if (nonOrthoCoeff < maxNonOrthoRatio)
{
orthogonal_ = true;
deleteDemandDrivenData(correctionVectors_);
}
else
{
orthogonal_ = false;
nonOrthCorrectionVectorsPtr_.reset(nullptr);
}
orthogonal_ = !bool(nonOrthCorrectionVectorsPtr_);
DebugInFunction
<< "Finished constructing non-orthogonal correction vectors"
<< endl;
@ -547,21 +581,26 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
<< "Constructing skew correction vectors"
<< endl;
skewCorrectionVectors_ = new edgeVectorField
const polyMesh& pMesh = mesh().mesh();
skewCorrectionVectorsPtr_.reset
(
IOobject
new edgeVectorField
(
"skewCorrectionVectors",
mesh()().pointsInstance(),
mesh()(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh(),
dimless
IOobject
(
"skewCorrectionVectors",
pMesh.pointsInstance(),
pMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh(),
dimensionedVector(dimless, Zero)
)
);
edgeVectorField& SkewCorrVecs = *skewCorrectionVectors_;
edgeVectorField& skewCorrVecs = *skewCorrectionVectorsPtr_;
// Set local references to mesh data
const areaVectorField& C = mesh().areaCentres();
@ -574,116 +613,125 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
const edgeList& edges = mesh().edges();
forAll(neighbour, edgeI)
forAll(neighbour, edgei)
{
const vector& P = C[owner[edgeI]];
const vector& N = C[neighbour[edgeI]];
const vector& S = points[edges[edgeI].start()];
vector e = edges[edgeI].vec(points);
const vector& P = C[owner[edgei]];
const vector& N = C[neighbour[edgei]];
const vector& S = points[edges[edgei].start()];
const scalar beta = magSqr((N - P)^e);
// (T:Eq. 5.4)
const vector d(N - P);
const vector e(edges[edgei].vec(points));
const vector de(d^e);
const scalar alpha = magSqr(de);
if (beta < ROOTVSMALL)
if (alpha < SMALL)
{
// Too small - skew correction remains zero
continue;
}
const scalar beta = -((d^(S - P)) & de)/alpha;
const scalar alpha = -(((N - P)^(S - P))&((N - P)^e))/beta;
// (T:Eq. 5.3)
const vector E(S + beta*e);
vector E(S + alpha*e);
SkewCorrVecs[edgeI] = Ce[edgeI] - E;
skewCorrVecs[edgei] = Ce[edgei] - E;
}
edgeVectorField::Boundary& bSkewCorrVecs =
SkewCorrVecs.boundaryFieldRef();
skewCorrVecs.boundaryFieldRef();
forAll(SkewCorrVecs.boundaryField(), patchI)
forAll(skewCorrVecs.boundaryField(), patchi)
{
faePatchVectorField& patchSkewCorrVecs = bSkewCorrVecs[patchI];
faePatchVectorField& patchSkewCorrVecs = bSkewCorrVecs[patchi];
if (patchSkewCorrVecs.coupled())
{
const labelUList& edgeFaces =
mesh().boundary()[patchI].edgeFaces();
mesh().boundary()[patchi].edgeFaces();
const edgeList::subList patchEdges =
mesh().boundary()[patchI].patchSlice(edges);
mesh().boundary()[patchi].patchSlice(edges);
vectorField ngbC(C.boundaryField()[patchI].patchNeighbourField());
vectorField ngbC(C.boundaryField()[patchi].patchNeighbourField());
forAll(patchSkewCorrVecs, edgeI)
forAll(patchSkewCorrVecs, edgei)
{
const vector& P = C[edgeFaces[edgeI]];
const vector& N = ngbC[edgeI];
const vector& S = points[patchEdges[edgeI].start()];
vector e = patchEdges[edgeI].vec(points);
const vector& P = C[edgeFaces[edgei]];
const vector& N = ngbC[edgei];
const vector& S = points[patchEdges[edgei].start()];
const scalar beta = magSqr((N - P)^e);
// (T:Eq. 5.4)
const vector d(N - P);
const vector e(patchEdges[edgei].vec(points));
const vector de(d^e);
const scalar alpha = magSqr(de);
if (beta < ROOTVSMALL)
if (alpha < SMALL)
{
// Too small - skew correction remains zero
continue;
}
const scalar beta = -((d^(S - P)) & de)/alpha;
const scalar alpha = -(((N - P)^(S - P))&((N - P)^e))/beta;
const vector E(S + beta*e);
vector E(S + alpha*e);
patchSkewCorrVecs[edgeI] =
Ce.boundaryField()[patchI][edgeI] - E;
patchSkewCorrVecs[edgei] =
Ce.boundaryField()[patchi][edgei] - E;
}
}
else
{
patchSkewCorrVecs = vector::zero;
}
}
#ifdef FA_SKEW_CORRECTION
scalar skewCoeff = 0.0;
constexpr scalar maxSkewRatio = 0.1;
scalar skewCoeff = 0;
// Calculating PN arc length
scalarField lPN(owner.size());
forAll(owner, edgeI)
forAll(own, edgei)
{
lPN[edgeI] =
const scalar magSkew = mag(skewCorrVecs[edgei]);
const scalar lPN =
mag
(
Ce[edgeI]
- SkewCorrVecs[edgeI]
- C[owner[edgeI]]
Ce[edgei]
- skewCorrVecs[edgei]
- C[owner[edgei]]
)
+ mag
(
C[neighbour[edgeI]]
- Ce[edgeI]
+ SkewCorrVecs[edgeI]
C[neighbour[edgei]]
- Ce[edgei]
+ skewCorrVecs[edgei]
);
}
if (lPN.size() > 0)
{
skewCoeff = max(mag(SkewCorrVecs.internalField())/mag(lPN));
const scalar ratio = magSkew/lPN;
if (skewCoeff < ratio)
{
skewCoeff = ratio;
if (skewCoeff > maxSkewRatio)
{
break;
}
}
}
DebugInFunction
<< "skew coefficient = " << skewCoeff << endl;
if (skewCoeff < 0.1)
if (skewCoeff < maxSkewRatio)
{
skew_ = false;
deleteDemandDrivenData(skewCorrectionVectors_);
}
else
{
skew_ = true;
skewCorrectionVectorsPtr_.reset(nullptr);
}
#endif
skew_ = bool(skewCorrectionVectorsPtr_);
DebugInFunction
<< "Finished constructing skew correction vectors"
<< endl;

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.
@ -70,19 +71,19 @@ class edgeInterpolation
// Demand-driven data
//- Geodesic distance between centroids of neighbour finite areas
mutable edgeScalarField* lPN_;
mutable unique_ptr<edgeScalarField> lPNptr_;
//- Central-differencing weighting factors
mutable edgeScalarField* weightingFactors_;
mutable unique_ptr<edgeScalarField> weightsPtr_;
//- Face-gradient difference factors
mutable edgeScalarField* differenceFactors_;
mutable unique_ptr<edgeScalarField> deltaCoeffsPtr_;
//- Non-orthogonality correction vectors
mutable edgeVectorField* correctionVectors_;
mutable unique_ptr<edgeVectorField> nonOrthCorrectionVectorsPtr_;
//- Skew correction vectors
mutable edgeVectorField* skewCorrectionVectors_;
mutable unique_ptr<edgeVectorField> skewCorrectionVectorsPtr_;
//- Is mesh orthogonal
mutable bool orthogonal_;
@ -93,6 +94,9 @@ class edgeInterpolation
// Private Member Functions
//- Return skewness correction per edge
const vector& skewCorr(const label edgeI) const;
//- Construct geodesic distance between P and N
void makeLPN() const;
@ -116,7 +120,7 @@ protected:
// Storage Management
//- Clear all geometry and addressing
void clearOut();
void clearOut() const;
public:
@ -138,10 +142,7 @@ public:
// Member Functions
//- Return mesh reference
const faMesh& mesh() const noexcept
{
return faMesh_;
}
const faMesh& mesh() const noexcept { return faMesh_; }
//- Return reference to PN geodesic distance
const edgeScalarField& lPN() const;
@ -153,7 +154,7 @@ public:
const edgeScalarField& deltaCoeffs() const;
//- Return reference to non-orthogonality correction vectors array
const edgeVectorField& correctionVectors() const;
const edgeVectorField& nonOrthCorrectionVector() const;
//- Return reference to skew vectors array
const edgeVectorField& skewCorrectionVectors() const;
@ -174,7 +175,7 @@ public:
// Storage Management
//- True if weights exist
bool hasWeights() const noexcept { return bool(weightingFactors_); }
bool hasWeights() const noexcept { return bool(weightsPtr_); }
};

View File

@ -163,19 +163,16 @@ Foam::edgeInterpolationScheme<Type>::interpolate
const labelUList& P = mesh.owner();
const labelUList& N = mesh.neighbour();
tmp<GeometricField<Type, faePatchField, edgeMesh>> tsf
auto tsf = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
(
new GeometricField<Type, faePatchField, edgeMesh>
IOobject
(
IOobject
(
"interpolate("+vf.name()+')',
vf.instance(),
vf.db()
),
mesh,
vf.dimensions()
)
"interpolate("+vf.name()+')',
vf.instance(),
vf.db()
),
mesh,
vf.dimensions()
);
GeometricField<Type, faePatchField, edgeMesh>& sf = tsf.ref();
@ -208,8 +205,8 @@ Foam::edgeInterpolationScheme<Type>::interpolate
if (vf.boundaryField()[pi].coupled())
{
label size = vf.boundaryField()[pi].patch().size();
label start = vf.boundaryField()[pi].patch().start();
const label size = vf.boundaryField()[pi].patch().size();
const label start = vf.boundaryField()[pi].patch().start();
Field<Type> pOwnVf = vf.boundaryField()[pi].patchInternalField();
Field<Type> pNgbVf = vf.boundaryField()[pi].patchNeighbourField();
@ -233,10 +230,6 @@ Foam::edgeInterpolationScheme<Type>::interpolate
+ pY[i]*transform(TN, pNgbVf[i])
);
}
// sf.boundaryFieldRef()[pi] =
// pLambda*vf.boundaryField()[pi].patchInternalField()
// + pY*vf.boundaryField()[pi].patchNeighbourField();
}
else
{
@ -279,19 +272,16 @@ Foam::edgeInterpolationScheme<Type>::interpolate
const labelUList& P = mesh.owner();
const labelUList& N = mesh.neighbour();
tmp<GeometricField<Type, faePatchField, edgeMesh>> tsf
auto tsf = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
(
new GeometricField<Type, faePatchField, edgeMesh>
IOobject
(
IOobject
(
"interpolate("+vf.name()+')',
vf.instance(),
vf.db()
),
mesh,
vf.dimensions()
)
"interpolate("+vf.name()+')',
vf.instance(),
vf.db()
),
mesh,
vf.dimensions()
);
GeometricField<Type, faePatchField, edgeMesh>& sf = tsf.ref();
@ -325,8 +315,8 @@ Foam::edgeInterpolationScheme<Type>::interpolate
if (vf.boundaryField()[pi].coupled())
{
label size = vfb[pi].patch().size();
label start = vfb[pi].patch().start();
const label size = vfb[pi].patch().size();
const label start = vfb[pi].patch().start();
Field<Type> pOwnVf(vfb[pi].patchInternalField());
Field<Type> pNgbVf(vfb[pi].patchNeighbourField());
@ -350,10 +340,6 @@ Foam::edgeInterpolationScheme<Type>::interpolate
+ (1 - pLambda[i])*transform(TN, pNgbVf[i])
);
}
// tsf().boundaryFieldRef()[pi] =
// pLambda*vf.boundaryField()[pi].patchInternalField()
// + (1 - pLambda)*vf.boundaryField()[pi].patchNeighbourField();
}
else
{
@ -395,19 +381,16 @@ Foam::edgeInterpolationScheme<Type>::euclidianInterpolate
const labelUList& P = mesh.owner();
const labelUList& N = mesh.neighbour();
tmp<GeometricField<Type, faePatchField, edgeMesh>> tsf
auto tsf = tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
(
new GeometricField<Type, faePatchField, edgeMesh>
IOobject
(
IOobject
(
"interpolate("+vf.name()+')',
vf.instance(),
vf.db()
),
mesh,
vf.dimensions()
)
"interpolate("+vf.name()+')',
vf.instance(),
vf.db()
),
mesh,
vf.dimensions()
);
GeometricField<Type, faePatchField, edgeMesh>& sf = tsf.ref();

View File

@ -48,6 +48,7 @@ SourceFiles
namespace Foam
{
// Forward Declarations
class faMesh;
/*---------------------------------------------------------------------------*\
@ -59,9 +60,9 @@ class edgeInterpolationScheme
:
public refCount
{
// Private data
// Private Data
//- Hold reference to mesh
//- Reference to the mesh
const faMesh& mesh_;
@ -136,10 +137,7 @@ public:
// Member Functions
//- Return mesh reference
const faMesh& mesh() const
{
return mesh_;
}
const faMesh& mesh() const noexcept { return mesh_; }
//- Return the face-interpolate of the given cell field

View File

@ -54,7 +54,11 @@ namespace Foam
class GammaWeight
{
scalar k_;
// Private Data
//- Model coefficient [0,1]
scalar k_;
public:
@ -72,7 +76,7 @@ public:
// Rescale k_ to be >= 0 and <= 0.5 (TVD conformant)
// and avoid the /0 when k_ = 0
k_ = max(k_/2.0, SMALL);
k_ = max(0.5*k_, SMALL);
}
@ -87,10 +91,7 @@ public:
const vector& d
) const
{
scalar magd = mag(d);
vector dHat = d/mag(d);
scalar gradf = (phiN - phiP)/magd;
const vector dHat(normalised(d));
scalar gradcf;
scalar udWeight;
@ -109,8 +110,10 @@ public:
// Stabilise for division
gradcf = stabilise(gradcf, SMALL);
scalar phict = 1 - 0.5*gradf/gradcf;
scalar limiter = min(max(phict/k_, 0), 1);
const scalar gradf = (phiN - phiP)/mag(d);
const scalar phict = 1 - 0.5*gradf/gradcf;
const scalar limiter = min(max(phict/k_, 0), 1);
return limiter*cdWeight + (1 - limiter)*udWeight;
}

View File

@ -62,9 +62,9 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
{
const faMesh& mesh = this->mesh();
tmp<edgeScalarField> tWeightingFactors
auto tWeightingFactors = tmp<edgeScalarField>::New
(
new edgeScalarField(mesh.edgeInterpolation::weights())
mesh.edgeInterpolation::weights()
);
edgeScalarField& weightingFactors = tWeightingFactors.ref();
@ -73,7 +73,8 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
tmp<areaScalarField> tvf = limiter(phi);
const areaScalarField& vf = tvf();
const areaVectorField gradc(fac::grad(vf));
tmp<areaVectorField> tgradc = fac::grad(vf);
const areaVectorField& gradc = tgradc.cref();
// edgeVectorField d
// (
@ -84,7 +85,7 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
// if (!mesh.orthogonal())
// {
// d -=
// mesh.edgeInterpolation::correctionVectors()
// mesh.edgeInterpolation::nonOrthCorrectionVector()
// /mesh.edgeInterpolation::deltaCoeffs();
// }
@ -108,6 +109,12 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
d.normalise();
d *= mesh.edgeInterpolation::lPN().internalField()[edge];
// Do not allow any mag(val) < SMALL
if (mag(d) < SMALL)
{
d = vector::uniform(SMALL);
}
weights[edge] =
this->weight
(
@ -124,80 +131,86 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
auto& bWeights = weightingFactors.boundaryFieldRef();
forAll(bWeights, patchI)
forAll(bWeights, patchi)
{
if (bWeights[patchI].coupled())
if (bWeights[patchi].coupled())
{
scalarField& pWeights = bWeights[patchI];
scalarField& pWeights = bWeights[patchi];
const scalarField& pEdgeFlux = edgeFlux_.boundaryField()[patchI];
const scalarField& pEdgeFlux = edgeFlux_.boundaryField()[patchi];
scalarField pVfP(vf.boundaryField()[patchI].patchInternalField());
scalarField pVfP(vf.boundaryField()[patchi].patchInternalField());
scalarField pVfN(vf.boundaryField()[patchI].patchNeighbourField());
scalarField pVfN(vf.boundaryField()[patchi].patchNeighbourField());
vectorField pGradcP
(
gradc.boundaryField()[patchI].patchInternalField()
gradc.boundaryField()[patchi].patchInternalField()
);
vectorField pGradcN
(
gradc.boundaryField()[patchI].patchNeighbourField()
gradc.boundaryField()[patchi].patchNeighbourField()
);
vectorField CP
(
mesh.areaCentres().boundaryField()[patchI].patchInternalField()
mesh.areaCentres().boundaryField()[patchi].patchInternalField()
);
vectorField CN
(
mesh.areaCentres().boundaryField()[patchI]
mesh.areaCentres().boundaryField()[patchi]
.patchNeighbourField()
);
vectorField nP
(
mesh.faceAreaNormals().boundaryField()[patchI]
mesh.faceAreaNormals().boundaryField()[patchi]
.patchInternalField()
);
vectorField nN
(
mesh.faceAreaNormals().boundaryField()[patchI]
mesh.faceAreaNormals().boundaryField()[patchi]
.patchNeighbourField()
);
scalarField pLPN
(
mesh.edgeInterpolation::lPN().boundaryField()[patchI]
mesh.edgeInterpolation::lPN().boundaryField()[patchi]
);
forAll(pWeights, edgeI)
forAll(pWeights, edgei)
{
vector d(CN[edgeI] - CP[edgeI]);
vector d(CN[edgei] - CP[edgei]);
if (pEdgeFlux[edgeI] > 0)
if (pEdgeFlux[edgei] > 0)
{
d.removeCollinear(nP[edgeI]);
d.removeCollinear(nP[edgei]);
}
else
{
d.removeCollinear(nN[edgeI]);
d.removeCollinear(nN[edgei]);
}
d.normalise();
d *= pLPN[edgeI];
d *= pLPN[edgei];
pWeights[edgeI] =
// Do not allow any mag(val) < SMALL
if (mag(d) < SMALL)
{
d = vector::uniform(SMALL);
}
pWeights[edgei] =
this->weight
(
pWeights[edgeI],
pEdgeFlux[edgeI],
pVfP[edgeI],
pVfN[edgeI],
pGradcP[edgeI],
pGradcN[edgeI],
pWeights[edgei],
pEdgeFlux[edgei],
pVfP[edgei],
pVfN[edgei],
pGradcP[edgei],
pGradcN[edgei],
d
);
}

View File

@ -56,8 +56,9 @@ class blendedEdgeInterpolation
public linearEdgeInterpolation<Type>,
public upwindEdgeInterpolation<Type>
{
// Private data
// Private Data
//- Blending factor
const scalar blendingFactor_;

View File

@ -60,6 +60,7 @@ class skewCorrectedEdgeInterpolation
{
// Private Data
//- Edge-interpolation scheme
tmp<edgeInterpolationScheme<Type>> tScheme_;
@ -132,9 +133,8 @@ public:
const edgeVectorField& scv = mesh.skewCorrectionVectors();
tmp<GeometricField<Type, faePatchField, edgeMesh>> tsfCorr
(
new GeometricField<Type, faePatchField, edgeMesh>
auto tsfCorr =
tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
(
IOobject
(
@ -144,9 +144,7 @@ public:
),
mesh,
dimensioned<Type>(vf.dimensions(), Zero)
)
);
);
GeometricField<Type, faePatchField, edgeMesh>& corr = tsfCorr.ref();
for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)

View File

@ -55,8 +55,9 @@ class upwindEdgeInterpolation
:
virtual public edgeInterpolationScheme<Type>
{
// Private data
// Private Data
//- Face flux
const edgeScalarField& faceFlux_;

View File

@ -1,13 +1,13 @@
/* Region models */
regionFaModel/regionFaModel.C
thermalShellModel/thermalShellModel.C
thermalShellModel/thermalShellModelNew.C
vibrationShellModel/vibrationShellModel.C
vibrationShellModel/vibrationShellModelNew.C
thermalShellModels/thermalShellModel/thermalShellModel.C
thermalShellModels/thermalShellModel/thermalShellModelNew.C
vibrationShellModels/vibrationShellModel/vibrationShellModel.C
vibrationShellModels/vibrationShellModel/vibrationShellModelNew.C
/* Shell models */
thermalShell/thermalShell.C
KirchhoffShell/KirchhoffShell.C
thermalShellModels/thermalShell/thermalShell.C
vibrationShellModels/KirchhoffShell/KirchhoffShell.C
/* Boundary conditions */
derivedFvPatchFields/thermalShell/thermalShellFvPatchScalarField.C
@ -15,31 +15,32 @@ derivedFvPatchFields/vibrationShell/vibrationShellFvPatchScalarField.C
/* Sub-Model */
liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
kinematic = liquidFilmModels/subModels/kinematic
$(kinematic)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
$(kinematic)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
$(kinematic)/filmTurbulenceModel/laminar/laminar.C
liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C
liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C
liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C
$(kinematic)/injectionModel/injectionModelList/injectionModelList.C
$(kinematic)/injectionModel/injectionModel/injectionModel.C
$(kinematic)/injectionModel/injectionModel/injectionModelNew.C
liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C
$(kinematic)/injectionModel/curvatureSeparation/curvatureSeparation.C
$(kinematic)/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C
liquidFilm/subModels/kinematic/force/forceList/forceList.C
liquidFilm/subModels/kinematic/force/force/force.C
liquidFilm/subModels/kinematic/force/force/forceNew.C
liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C
liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C
$(kinematic)/force/forceList/forceList.C
$(kinematic)/force/force/force.C
$(kinematic)/force/force/forceNew.C
$(kinematic)/force/contactAngleForces/contactAngleForce/contactAngleForce.C
$(kinematic)/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C
liquidFilm/subModels/filmSubModelBase.C
liquidFilmModels/subModels/filmSubModelBase.C
liquidFilm/liquidFilmBase.C
liquidFilm/liquidFilmBaseNew.C
liquidFilmModels/liquidFilmBase.C
liquidFilmModels/liquidFilmBaseNew.C
liquidFilm/liquidFilmModel/liquidFilmModel.C
liquidFilm/kinematicThinFilm/kinematicThinFilm.C
derivedFvPatchFields/filmShell/velocityFilmShellFvPatchVectorField.C
liquidFilmModels/liquidFilmModel/liquidFilmModel.C
liquidFilmModels/kinematicThinFilm/kinematicThinFilm.C
derivedFvPatchFields/velocityFilmShell/velocityFilmShellFvPatchVectorField.C
functionObjects/setTimeStep/setTimeStepFaRegionsFunctionObject.C

View File

@ -44,8 +44,8 @@ thermalShellFvPatchScalarField::thermalShellFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchField<scalar>(p, iF),
baffle_(nullptr),
fixedValueFvPatchScalarField(p, iF),
bafflePtr_(nullptr),
dict_()
{}
@ -58,14 +58,14 @@ thermalShellFvPatchScalarField::thermalShellFvPatchScalarField
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchField<scalar>
fixedValueFvPatchScalarField
(
ptf,
p,
iF,
mapper
),
baffle_(nullptr),
bafflePtr_(nullptr),
dict_(ptf.dict_)
{}
@ -77,8 +77,8 @@ thermalShellFvPatchScalarField::thermalShellFvPatchScalarField
const dictionary& dict
)
:
fixedValueFvPatchField<scalar>(p, iF, dict),
baffle_(nullptr),
fixedValueFvPatchScalarField(p, iF, dict),
bafflePtr_(nullptr),
dict_
(
// Copy dictionary, but without "heavy" data chunks
@ -94,9 +94,9 @@ thermalShellFvPatchScalarField::thermalShellFvPatchScalarField
)
)
{
if (!baffle_)
if (!bafflePtr_)
{
baffle_.reset(baffleType::New(p.boundaryMesh().mesh(), dict_));
bafflePtr_.reset(baffleType::New(p.boundaryMesh().mesh(), dict_));
}
}
@ -107,8 +107,8 @@ thermalShellFvPatchScalarField::thermalShellFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchField<scalar>(ptf, iF),
baffle_(nullptr),
fixedValueFvPatchScalarField(ptf, iF),
bafflePtr_(nullptr),
dict_(ptf.dict_)
{}
@ -122,19 +122,19 @@ void thermalShellFvPatchScalarField::updateCoeffs()
return;
}
baffle_->evolve();
bafflePtr_->evolve();
scalarField& pfld = *this;
baffle_->vsm().mapToVolumePatch(baffle_->T(), pfld, patch().index());
bafflePtr_->vsm().mapToVolumePatch(bafflePtr_->T(), pfld, patch().index());
fixedValueFvPatchField<scalar>::updateCoeffs();
fixedValueFvPatchScalarField::updateCoeffs();
}
void thermalShellFvPatchScalarField::write(Ostream& os) const
{
fixedValueFvPatchField<scalar>::write(os);
fixedValueFvPatchScalarField::write(os);
dict_.write(os, false);
}

View File

@ -41,17 +41,17 @@ Usage
\verbatim
<masterPatchName>
{
// Mandatory entries (unmodifiable)
// Mandatory entries
type compressible::thermalShell;
// Mandatory/Optional (inherited) entries
// Inherited entries
...
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Dflt
Property | Description | Type | Reqd | Deflt
type | Type name: compressible::thermalShell | word | yes | -
\endtable
@ -88,7 +88,7 @@ namespace compressible
class thermalShellFvPatchScalarField
:
public fixedValueFvPatchField<scalar>
public fixedValueFvPatchScalarField
{
// Typedefs
@ -99,7 +99,7 @@ class thermalShellFvPatchScalarField
// Private Data
//- The thermal baffle
autoPtr<baffleType> baffle_;
autoPtr<baffleType> bafflePtr_;
//- Dictionary
dictionary dict_;

View File

@ -42,8 +42,8 @@ velocityFilmShellFvPatchVectorField::velocityFilmShellFvPatchVectorField
const DimensionedField<vector, volMesh>& iF
)
:
mixedFvPatchField<vector>(p, iF),
baffle_(nullptr),
mixedFvPatchVectorField(p, iF),
bafflePtr_(nullptr),
dict_(),
curTimeIndex_(-1),
zeroWallVelocity_(true)
@ -62,14 +62,14 @@ velocityFilmShellFvPatchVectorField::velocityFilmShellFvPatchVectorField
const fvPatchFieldMapper& mapper
)
:
mixedFvPatchField<vector>
mixedFvPatchVectorField
(
ptf,
p,
iF,
mapper
),
baffle_(nullptr),
bafflePtr_(nullptr),
dict_(ptf.dict_),
curTimeIndex_(-1),
zeroWallVelocity_(true)
@ -83,8 +83,8 @@ velocityFilmShellFvPatchVectorField::velocityFilmShellFvPatchVectorField
const dictionary& dict
)
:
mixedFvPatchField<vector>(p, iF),
baffle_(nullptr),
mixedFvPatchVectorField(p, iF),
bafflePtr_(nullptr),
dict_
(
// Copy dictionary, but without "heavy" data chunks
@ -119,9 +119,9 @@ velocityFilmShellFvPatchVectorField::velocityFilmShellFvPatchVectorField
valueFraction() = 1;
}
if (!baffle_)
if (!bafflePtr_)
{
baffle_.reset(baffleType::New(p.boundaryMesh().mesh(), dict_));
bafflePtr_.reset(baffleType::New(p.boundaryMesh().mesh(), dict_));
}
}
@ -132,8 +132,8 @@ velocityFilmShellFvPatchVectorField::velocityFilmShellFvPatchVectorField
const DimensionedField<vector, volMesh>& iF
)
:
mixedFvPatchField<vector>(ptf, iF),
baffle_(nullptr),
mixedFvPatchVectorField(ptf, iF),
bafflePtr_(nullptr),
dict_(ptf.dict_),
curTimeIndex_(-1),
zeroWallVelocity_(true)
@ -142,7 +142,6 @@ velocityFilmShellFvPatchVectorField::velocityFilmShellFvPatchVectorField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void velocityFilmShellFvPatchVectorField::updateCoeffs()
{
if (this->updated())
@ -150,14 +149,19 @@ void velocityFilmShellFvPatchVectorField::updateCoeffs()
return;
}
// Execute the change to the openFraction only once per time-step
// Execute the change only once per time-step
if (curTimeIndex_ != this->db().time().timeIndex())
{
baffle_->evolve();
bafflePtr_->evolve();
vectorField& pfld = *this;
baffle_->vsm().mapToVolumePatch(baffle_->Us(), pfld, patch().index());
bafflePtr_->vsm().mapToVolumePatch
(
bafflePtr_->Us(),
pfld,
patch().index()
);
refGrad() = Zero;
valueFraction() = 1;
@ -174,7 +178,7 @@ void velocityFilmShellFvPatchVectorField::updateCoeffs()
curTimeIndex_ = this->db().time().timeIndex();
}
mixedFvPatchField<vector>::updateCoeffs();
mixedFvPatchVectorField::updateCoeffs();
}

View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -30,25 +30,30 @@ Group
grpLiquidFilmBoundaryConditions
Description
This boundary condition provides a coupled velocity condition between
a primary region (3D mesh) and a liquid-film model (2D mesh).
Usage
Example of the boundary condition specification:
\verbatim
<patchName>
{
// Mandatory entries
type velocityFilmShell;
liquidFilmModel <word>;
active true;
infoOutput true;
// Optional entries
zeroWallVelocity <bool>;
U U;
pRef 1e5;
T0 300;
deltaWet 1e-4;
h0 1e-8;
zeroWallVelocity true;
// Inherited entries
...
active <bool>;
infoOutput <bool>;
U <word>;
pRef <scalar>;
T0 <scalar>;
deltaWet <scalar>;
h0 <scalar>;
thermo
{
@ -56,7 +61,6 @@ Usage
}
turbulence laminar;
laminarCoeffs
{
friction ManningStrickler; // Wall friction model
@ -77,8 +81,6 @@ Usage
}
region film;
liquidFilmModel kinematicThinFilm;
value uniform (0 0 0);
}
\endverbatim
@ -87,28 +89,32 @@ Usage
\table
Property | Description | Type | Reqd | Deflt
type | Type name: velocityFilmShell | word | yes | -
liquidFilmModel | Film model | word | yes | -
zeroWallVelocity | Flag to fix zero U for primary flow <!--
--> | bool | no | true
U | Name of the primary U | word | yes | -
pRef | Reference pressure for thermo | scalar | yes | -
T0 | Film initial temperature | scalar | no | READ
thermo | Flow thermo | wordRes | yes | -
zeroWallVelocity | Flag to fix zero U for primary flow <!--
--> | bool | no | true
turbulence | Type of film turbulence model | word | yes | -
injectionModels | Lagrangian injection | | no | -
forces | Film force models | wordRes | no | -
deltaWet | Wet film thickness | scalar | no | 1e-4
h0 | Numerical minimum thickness | scalar | no | 1e-7
region | Name of the 2D region | word | yes | -
liquidFilmModel | Film model | word | yes | -
\endtable
The inherited entries are elaborated in:
- \link liquidFilmBase.H \endlink
- \link mixedFvPatchFields.H \endlink
SourceFiles
velocityFilmShellFvPatchVectorField.C
\*---------------------------------------------------------------------------*/
#ifndef velocityFilmShellFvPatchVectorField_H
#define velocityFilmShellFvPatchVectorField_H
#ifndef Foam_velocityFilmShellFvPatchVectorField_H
#define Foam_velocityFilmShellFvPatchVectorField_H
#include "autoPtr.H"
#include "liquidFilmBase.H"
@ -125,7 +131,7 @@ namespace Foam
class velocityFilmShellFvPatchVectorField
:
public mixedFvPatchField<vector>
public mixedFvPatchVectorField
{
// Typedefs
@ -135,8 +141,8 @@ class velocityFilmShellFvPatchVectorField
// Private Data
//- The liquid film
autoPtr<baffleType> baffle_;
//- The liquid film model
autoPtr<baffleType> bafflePtr_;
//- Dictionary
mutable dictionary dict_;
@ -144,7 +150,7 @@ class velocityFilmShellFvPatchVectorField
//- Time index to evolve the film
label curTimeIndex_;
//- Zero wall velocity. Fix U to zero or to film U
//- Flag to set velocity to zero or film velocity
bool zeroWallVelocity_;

View File

@ -42,8 +42,8 @@ vibrationShellFvPatchScalarField::vibrationShellFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
mixedFvPatchField<scalar>(p, iF),
baffle_(nullptr),
mixedFvPatchScalarField(p, iF),
bafflePtr_(nullptr),
dict_()
{
refValue() = Zero;
@ -60,14 +60,14 @@ vibrationShellFvPatchScalarField::vibrationShellFvPatchScalarField
const fvPatchFieldMapper& mapper
)
:
mixedFvPatchField<scalar>
mixedFvPatchScalarField
(
ptf,
p,
iF,
mapper
),
baffle_(nullptr),
bafflePtr_(nullptr),
dict_(ptf.dict_)
{}
@ -79,8 +79,8 @@ vibrationShellFvPatchScalarField::vibrationShellFvPatchScalarField
const dictionary& dict
)
:
mixedFvPatchField<scalar>(p, iF),
baffle_(nullptr),
mixedFvPatchScalarField(p, iF),
bafflePtr_(nullptr),
dict_
(
// Copy dictionary, but without "heavy" data chunks
@ -113,9 +113,9 @@ vibrationShellFvPatchScalarField::vibrationShellFvPatchScalarField
valueFraction() = 1;
}
if (!baffle_)
if (!bafflePtr_)
{
baffle_.reset(baffleType::New(p.boundaryMesh().mesh(), dict_));
bafflePtr_.reset(baffleType::New(p.boundaryMesh().mesh(), dict_));
}
}
@ -126,8 +126,8 @@ vibrationShellFvPatchScalarField::vibrationShellFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
mixedFvPatchField<scalar>(ptf, iF),
baffle_(nullptr),
mixedFvPatchScalarField(ptf, iF),
bafflePtr_(nullptr),
dict_(ptf.dict_)
{}
@ -144,28 +144,33 @@ void vibrationShellFvPatchScalarField::updateCoeffs()
const auto& transportProperties =
db().lookupObject<IOdictionary>("transportProperties");
dimensionedScalar rho("rho", dimDensity, transportProperties);
const dimensionedScalar rho("rho", dimDensity, transportProperties);
baffle_->evolve();
bafflePtr_->evolve();
// rho * acceleration
refGrad() = Zero; // safety (for any unmapped values)
baffle_->vsm().mapToVolumePatch(baffle_->a(), refGrad(), patch().index());
bafflePtr_->vsm().mapToVolumePatch
(
bafflePtr_->a(),
refGrad(),
patch().index()
);
refGrad() *= rho.value();
refValue() = Zero;
valueFraction() = Zero;
mixedFvPatchField<scalar>::updateCoeffs();
mixedFvPatchScalarField::updateCoeffs();
}
void vibrationShellFvPatchScalarField::write(Ostream& os) const
{
mixedFvPatchField<scalar>::write(os);
mixedFvPatchScalarField::write(os);
dict_.write(os, false);
}

View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019-2020 OpenCFD Ltd.
Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -30,23 +30,25 @@ Group
grpVibrationBoundaryConditions
Description
This boundary condition provides a coupled acoustic pressure condition
between a primary region (3D mesh) and a vibration shell model (2D mesh).
Usage
Example of the boundary condition specification:
\verbatim
<masterPatchName>
{
// Mandatory entries (unmodifiable)
// Mandatory entries
type vibrationShell;
// Mandatory/Optional (inherited) entries
// Inherited entries
...
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Dflt
Property | Description | Type | Reqd | Deflt
type | Type name: vibrationShell | word | yes | -
\endtable
@ -59,8 +61,8 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef vibrationShellFvPatchScalarField_H
#define vibrationShellFvPatchScalarField_H
#ifndef Foam_vibrationShellFvPatchScalarField_H
#define Foam_vibrationShellFvPatchScalarField_H
#include "autoPtr.H"
#include "vibrationShellModel.H"
@ -77,7 +79,7 @@ namespace Foam
class vibrationShellFvPatchScalarField
:
public mixedFvPatchField<scalar>
public mixedFvPatchScalarField
{
// Typedefs
@ -87,8 +89,8 @@ class vibrationShellFvPatchScalarField
// Private Data
//- The vibration shell
autoPtr<baffleType> baffle_;
//- The vibration shell model
autoPtr<baffleType> bafflePtr_;
//- Dictionary
dictionary dict_;

View File

@ -75,10 +75,7 @@ void kinematicThinFilm::preEvolveRegion()
void kinematicThinFilm::evolveRegion()
{
if (debug)
{
InfoInFunction << endl;
}
DebugInFunction << endl;
const areaVectorField& ns = regionMesh().faceAreaNormals();
@ -86,7 +83,7 @@ void kinematicThinFilm::evolveRegion()
phi2s_ = fac::interpolate(h_)*phif_;
for (int oCorr=1; oCorr<=nOuterCorr_; oCorr++)
for (int oCorr=1; oCorr<=nOuterCorr_; ++oCorr)
{
pf_.storePrevIter();
@ -96,7 +93,7 @@ void kinematicThinFilm::evolveRegion()
+ fam::div(phi2s_, Uf_)
==
gs*h_
+ turbulence_->Su(Uf_)
+ turbulencePtr_->Su(Uf_)
+ faOptions()(h_, Uf_, sqr(dimVelocity))
+ forces_.correct(Uf_)
+ USp_
@ -108,25 +105,27 @@ void kinematicThinFilm::evolveRegion()
if (momentumPredictor_)
{
solve(UsEqn == - fac::grad(pf_*h_)/rho_ + pf_*fac::grad(h_)/rho_);
solve(UsEqn == -fac::grad(pf_*h_)/rho_ + pf_*fac::grad(h_)/rho_);
}
for (int corr=1; corr<=nCorr_; corr++)
for (int corr=1; corr<=nCorr_; ++corr)
{
areaScalarField UsA(UsEqn.A());
tmp<areaScalarField> tUsA = UsEqn.A();
Uf_ = UsEqn.H()/UsA;
Uf_ = UsEqn.H()/tUsA();
Uf_.correctBoundaryConditions();
faOptions().correct(Uf_);
const areaScalarField invRhoUsA(scalar(1)/(rho_*tUsA));
phif_ =
(fac::interpolate(Uf_) & regionMesh().Le())
- fac::interpolate(1.0/(rho_*UsA))
- fac::interpolate(invRhoUsA)
* fac::lnGrad(pf_*h_)*regionMesh().magLe()
+ fac::interpolate(pf_/(rho_*UsA))
+ fac::interpolate(pf_*invRhoUsA)
* fac::lnGrad(h_)*regionMesh().magLe();
for (int nFilm=1; nFilm<=nFilmCorr_; nFilm++)
for (int nFilm=1; nFilm<=nFilmCorr_; ++nFilm)
{
faScalarMatrix hEqn
(
@ -155,17 +154,15 @@ void kinematicThinFilm::evolveRegion()
pf_.correctBoundaryConditions();
pf_.relax();
Uf_ -= (1.0/(rho_*UsA))*fac::grad(pf_*h_)
- (pf_/(rho_*UsA))*fac::grad(h_);
Uf_ -= invRhoUsA*fac::grad(pf_*h_)
- (pf_*invRhoUsA)*fac::grad(h_);
Uf_.correctBoundaryConditions();
faOptions().correct(Uf_);
}
}
Info<< "Film h min/max = " << min(h_).value() << ", "
<< max(h_).value() << endl;
Info<< "Film U min/max = " << max(mag(Uf_)).value() << endl;
Info<< "Film h min/max = " << gMinMax(h_) << nl
<< "Film mag(U) min/max = " << gMinMaxMag(Uf_) << endl;
}
@ -178,7 +175,7 @@ void kinematicThinFilm::postEvolveRegion()
correctThermoFields();
// Correct turbulence
turbulence_->correct();
turbulencePtr_->correct();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -29,14 +29,34 @@ Class
Description
Thin film model.
Usage
Example of the boundary condition specification:
\verbatim
{
// Mandatory entries
liquidFilmModel kinematicThinFilm;
// Inherited entries
...
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
liquidFilmModel | Type name: kinematicThinFilm | dict | yes | -
\endtable
The inherited entries are elaborated in:
- \link liquidFilmModel.H \endlink
SourceFiles
kinematicThinFilm.C
kinematicThinFilmI.H
\*---------------------------------------------------------------------------*/
#ifndef kinematicThinFilm_H
#define kinematicThinFilm_H
#ifndef areaSurfaceFilmModels_kinematicThinFilm_H
#define areaSurfaceFilmModels_kinematicThinFilm_H
#include "volFieldsFwd.H"
#include "liquidFilmModel.H"
@ -89,16 +109,16 @@ public:
// Member Functions
// Evolution
// Evolution
//- Pre-evolve film
virtual void preEvolveRegion();
//- Pre-evolve film
virtual void preEvolveRegion();
//- Evolve the film
virtual void evolveRegion();
//- Evolve the film
virtual void evolveRegion();
//- Post-evolve film
virtual void postEvolveRegion();
//- Post-evolve film
virtual void postEvolveRegion();
};

View File

@ -60,7 +60,6 @@ liquidFilmBase::liquidFilmBase
)
:
regionFaModel(mesh, liquidFilmName, modelType, dict, true),
momentumPredictor_
(
this->solution().subDict("PIMPLE").get<bool>("momentumPredictor")
@ -74,17 +73,11 @@ liquidFilmBase::liquidFilmBase
(
this->solution().subDict("PIMPLE").get<label>("nFilmCorr")
),
h0_("h0", dimLength, 1e-7, dict),
deltaWet_("deltaWet", dimLength, 1e-4, dict),
UName_(dict.get<word>("U")),
pName_(dict.lookupOrDefault<word>("p", word::null)),
pName_(dict.getOrDefault<word>("p", word::null)),
pRef_(dict.get<scalar>("pRef")),
h_
(
IOobject
@ -97,7 +90,6 @@ liquidFilmBase::liquidFilmBase
),
regionMesh()
),
Uf_
(
IOobject
@ -148,7 +140,6 @@ liquidFilmBase::liquidFilmBase
),
fac::interpolate(Uf_) & regionMesh().Le()
),
phi2s_
(
IOobject
@ -161,8 +152,6 @@ liquidFilmBase::liquidFilmBase
),
fac::interpolate(h_*Uf_) & regionMesh().Le()
),
gn_
(
IOobject
@ -176,9 +165,7 @@ liquidFilmBase::liquidFilmBase
regionMesh(),
dimensionedScalar(dimAcceleration, Zero)
),
g_(meshObjects::gravity::New(primaryMesh().time())),
massSource_
(
IOobject
@ -191,7 +178,6 @@ liquidFilmBase::liquidFilmBase
dimensionedScalar(dimMass, Zero),
calculatedFvPatchField<scalar>::typeName
),
momentumSource_
(
IOobject
@ -204,7 +190,6 @@ liquidFilmBase::liquidFilmBase
dimensionedVector(dimPressure, Zero),
calculatedFvPatchField<vector>::typeName
),
pnSource_
(
IOobject
@ -218,21 +203,7 @@ liquidFilmBase::liquidFilmBase
calculatedFvPatchField<scalar>::typeName
),
energySource_
(
IOobject
(
"energySource",
primaryMesh().time().timeName(),
primaryMesh()
),
primaryMesh(),
dimensionedScalar(dimEnergy, Zero),
calculatedFvPatchField<scalar>::typeName
),
addedMassTotal_(0),
faOptions_(Foam::fa::options::New(primaryMesh()))
{
const areaVectorField& ns = regionMesh().faceAreaNormals();
@ -256,18 +227,15 @@ liquidFilmBase::~liquidFilmBase()
scalar liquidFilmBase::CourantNumber() const
{
scalar CoNum = 0.0;
scalar velMag = 0.0;
edgeScalarField SfUfbyDelta
const edgeScalarField SfUfbyDelta
(
regionMesh().edgeInterpolation::deltaCoeffs()*mag(phif_)
);
CoNum =
scalar CoNum =
max(SfUfbyDelta/regionMesh().magLe()).value()*time().deltaT().value();
velMag = max(mag(phif_)/regionMesh().magLe()).value();
scalar velMag = max(mag(phif_)/regionMesh().magLe()).value();
reduce(CoNum, maxOp<scalar>());
reduce(velMag, maxOp<scalar>());
@ -281,19 +249,16 @@ scalar liquidFilmBase::CourantNumber() const
Foam::tmp<Foam::areaVectorField> liquidFilmBase::Uw() const
{
tmp<areaVectorField> tUw
auto tUw = tmp<areaVectorField>::New
(
new areaVectorField
IOobject
(
IOobject
(
"tUw",
primaryMesh().time().timeName(),
primaryMesh()
),
regionMesh(),
dimensionedVector(dimVelocity, Zero)
)
"tUw",
primaryMesh().time().timeName(),
primaryMesh()
),
regionMesh(),
dimensionedVector(dimVelocity, Zero)
);
auto& Uw = tUw.ref();
@ -336,19 +301,16 @@ Foam::tmp<Foam::areaVectorField> liquidFilmBase::Uw() const
Foam::tmp<Foam::areaVectorField> liquidFilmBase::Us() const
{
tmp<areaVectorField> tUs
auto tUs = tmp<areaVectorField>::New
(
new areaVectorField
IOobject
(
IOobject
(
"tUs",
primaryMesh().time().timeName(),
primaryMesh()
),
regionMesh(),
dimensionedVector(dimVelocity, Zero)
)
"tUs",
primaryMesh().time().timeName(),
primaryMesh()
),
regionMesh(),
dimensionedVector(dimVelocity, Zero)
);
// Uf quadratic profile
@ -360,22 +322,18 @@ Foam::tmp<Foam::areaVectorField> liquidFilmBase::Us() const
Foam::tmp<Foam::areaVectorField> liquidFilmBase::Up() const
{
const volVectorField& Uprimary =
primaryMesh().lookupObject<volVectorField>(UName_);
const auto& Uprimary = primaryMesh().lookupObject<volVectorField>(UName_);
tmp<areaVectorField> tUp
auto tUp = tmp<areaVectorField>::New
(
new areaVectorField
IOobject
(
IOobject
(
"tUp",
primaryMesh().time().timeName(),
primaryMesh()
),
regionMesh(),
dimensionedVector(dimVelocity, Zero)
)
"tUp",
primaryMesh().time().timeName(),
primaryMesh()
),
regionMesh(),
dimensionedVector(dimVelocity, Zero)
);
areaVectorField& Up = tUp.ref();
@ -411,19 +369,16 @@ Foam::tmp<Foam::areaVectorField> liquidFilmBase::Up() const
tmp<areaScalarField> liquidFilmBase::pg() const
{
tmp<areaScalarField> tpg
auto tpg = tmp<areaScalarField>::New
(
new areaScalarField
IOobject
(
IOobject
(
"tpg",
primaryMesh().time().timeName(),
primaryMesh()
),
regionMesh(),
dimensionedScalar(dimPressure, Zero)
)
"tpg",
primaryMesh().time().timeName(),
primaryMesh()
),
regionMesh(),
dimensionedScalar(dimPressure, Zero)
);
auto& pfg = tpg.ref();
@ -434,8 +389,6 @@ tmp<areaScalarField> liquidFilmBase::pg() const
primaryMesh().lookupObject<volScalarField>(pName_),
pfg.primitiveFieldRef()
);
//pfg -= pRef_;
}
return tpg;
@ -444,19 +397,16 @@ tmp<areaScalarField> liquidFilmBase::pg() const
tmp<areaScalarField> liquidFilmBase::alpha() const
{
tmp<areaScalarField> talpha
auto talpha = tmp<areaScalarField>::New
(
new areaScalarField
IOobject
(
IOobject
(
"talpha",
primaryMesh().time().timeName(),
primaryMesh()
),
regionMesh(),
dimensionedScalar(dimless, Zero)
)
"talpha",
primaryMesh().time().timeName(),
primaryMesh()
),
regionMesh(),
dimensionedScalar(dimless, Zero)
);
auto& alpha = talpha.ref();
@ -505,72 +455,6 @@ void liquidFilmBase::postEvolveRegion()
}
Foam::fa::options& liquidFilmBase::faOptions()
{
return faOptions_;
}
const areaVectorField& liquidFilmBase::Uf() const
{
return Uf_;
}
const areaScalarField& liquidFilmBase::gn() const
{
return gn_;
}
const uniformDimensionedVectorField& liquidFilmBase::g() const
{
return g_;
}
const areaScalarField& liquidFilmBase::h() const
{
return h_;
}
const edgeScalarField& liquidFilmBase::phif() const
{
return phif_;
}
const edgeScalarField& liquidFilmBase::phi2s() const
{
return phi2s_;
}
const dimensionedScalar& liquidFilmBase::h0() const
{
return h0_;
}
const regionFaModel& liquidFilmBase::region() const
{
return *this;
}
scalar liquidFilmBase::pRef()
{
return pRef_;
}
word liquidFilmBase::UName() const
{
return UName_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace areaSurfaceFilmModels

View File

@ -27,7 +27,44 @@ Class
Foam::regionModels::liquidFilmBase
Description
Base class for thermal 2D shells
Base class for liquid-film models.
Usage
Example of the boundary condition specification:
\verbatim
<patchName>
{
// Mandatory entries
U <word>;
pRef <scalar>;
// Optional entries
h0 <scalar>;
deltaWet <scalar>;
p <word>;
// Inherited entries
...
momentumPredictor <int>;
nOuterCorr <label>;
nCorr <label>;
nFilmCorr <label>;
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
U | Name of velocity field | word | yes | -
pRef | Reference absolute pressure | scalar | yes | -
h0 | Smallest film thickness | scalar | no | 1e-7
deltaWet | Film thickness beyond which face is assumed to be wet <!--
--> | scalar | no | 1e-4
p | Name of pressure field | word | no | null
\endtable
The inherited entries are elaborated in:
- \link regionFaModel.H \endlink
SourceFiles
liquidFilmBase.C
@ -55,7 +92,7 @@ namespace areaSurfaceFilmModels
{
/*---------------------------------------------------------------------------*\
Class liquidFilmBase Declaration
Class liquidFilmBase Declaration
\*---------------------------------------------------------------------------*/
class liquidFilmBase
@ -68,13 +105,13 @@ protected:
// Solution parameters
//- Momentum predictor
//- Flat to enable momentum predictor
Switch momentumPredictor_;
//- Number of outer correctors
label nOuterCorr_;
//- Number of PISO-like correctors
//- Number of PISO-like inner correctors
label nCorr_;
//- Number of film thickness correctors
@ -83,10 +120,10 @@ protected:
//- Cumulative continuity error
scalar cumulativeContErr_;
//- Smallest numerical thickness
//- Smallest film thickness
dimensionedScalar h0_;
//- Delta wet for sub-models
//- Film thickness beyond which face is assumed to be wet
dimensionedScalar deltaWet_;
@ -102,7 +139,7 @@ protected:
// Fields
//- Film hight
//- Film height
areaScalarField h_;
//- Film velocity
@ -138,9 +175,6 @@ protected:
//- Normal pressure by particles
volScalarField pnSource_;
//- Energy
volScalarField energySource_;
//- Total mass added
scalar addedMassTotal_;
@ -173,7 +207,6 @@ public:
// Constructors
//- Construct from type name and mesh and dict
liquidFilmBase
(
@ -209,115 +242,109 @@ public:
virtual scalar CourantNumber() const;
// Helper functions
// Helper functions
//- Wall velocity
tmp<areaVectorField> Uw() const;
//- Return wall velocity
tmp<areaVectorField> Uw() const;
//- Film surface film velocity
tmp<areaVectorField> Us() const;
//- Return film-surface velocity
tmp<areaVectorField> Us() const;
//- Primary region velocity at film hight. Assume the film to be
// in the viscous sub-layer
tmp<areaVectorField> Up() const;
//- Return primary region velocity at film height
// Assuming the film to be within the viscous sub-layer
tmp<areaVectorField> Up() const;
//- Map primary static pressure
tmp<areaScalarField> pg() const;
//- Return primary static pressure
tmp<areaScalarField> pg() const;
//- Wet indicator using h0
tmp<areaScalarField> alpha() const;
//- Return wet indicator using h0
tmp<areaScalarField> alpha() const;
// Access functions
// Access
//- Return faOptions
Foam::fa::options& faOptions();
//- Return faOptions
Foam::fa::options& faOptions() { return faOptions_; }
//- Access const reference Uf
const areaVectorField& Uf() const;
//- Return const reference to Uf
const areaVectorField& Uf() const noexcept { return Uf_; }
//- Access const reference gn
const areaScalarField& gn() const;
//- Return const reference to gn
const areaScalarField& gn() const noexcept { return gn_; }
//- Gravity
const uniformDimensionedVectorField& g() const;
//- Return gravity
const uniformDimensionedVectorField& g() const noexcept { return g_; }
//- Access const reference h
const areaScalarField& h() const;
//- Return const reference to h
const areaScalarField& h() const noexcept { return h_; }
//- Access to momentum flux
const edgeScalarField& phif() const;
//- Return const reference to phif
const edgeScalarField& phif() const noexcept { return phif_; }
//- Access continuity flux
const edgeScalarField& phi2s() const;
//- Return const reference to phi2s
const edgeScalarField& phi2s() const noexcept { return phi2s_; }
//- Return const reference to h0
const dimensionedScalar& h0() const noexcept { return h0_; }
//- Return const reference to this region
const regionFaModel& region() const noexcept { return *this; }
//- Return pRef
scalar pRef() const noexcept { return pRef_; }
//- Return UName
word UName() const noexcept { return UName_; }
//- Return h0
const dimensionedScalar& h0() const;
// Transfer fields - to the primary region (lagragian injection)
//- Access to this region
const regionFaModel& region() const;
//- Return the film mass available for transfer to cloud
virtual const volScalarField& cloudMassTrans() const = 0;
//- Access to pRef
scalar pRef();
//- Name of the U field
word UName() const;
//- Return the parcel diameters originating from film to cloud
virtual const volScalarField& cloudDiameterTrans() const = 0;
// Transfer fields - to the primary region (lagragian injection)
// Evolution
//- Return mass transfer source - Eulerian phase only
//virtual tmp<volScalarField> primaryMassTrans() const = 0;
//- Pre-evolve film
virtual void preEvolveRegion();
//- Return the film mass available for transfer to cloud
virtual const volScalarField& cloudMassTrans() const = 0;
//- Return the parcel diameters originating from film to cloud
virtual const volScalarField& cloudDiameterTrans() const = 0;
//- Post-evolve film
virtual void postEvolveRegion();
// Thermo variables
// Evolution
//- Access const reference mu
virtual const areaScalarField& mu() const = 0;
//- Pre-evolve film
virtual void preEvolveRegion();
//- Access const reference rho
virtual const areaScalarField& rho() const = 0;
//- Post-evolve film
virtual void postEvolveRegion();
//- Access const reference sigma
virtual const areaScalarField& sigma() const = 0;
//- Access const reference Tf
virtual const areaScalarField& Tf() const = 0;
//- Access const reference Cp
virtual const areaScalarField& Cp() const = 0;
// Thermo variables
// External hook to add sources and mass exchange variables
//- Access const reference mu
virtual const areaScalarField& mu() const = 0;
//- Access const reference rho
virtual const areaScalarField& rho() const = 0;
//- Access const reference sigma
virtual const areaScalarField& sigma() const = 0;
//- Access const reference Tf
virtual const areaScalarField& Tf() const = 0;
//- Access const reference Cp
virtual const areaScalarField& Cp() const = 0;
// External hook to add sources and mass exchange variables
//- Add sources
virtual void addSources
(
const label patchi, // patchi on primary region
const label facei, // facei of patchi
const scalar massSource, // [kg]
const vector& momentumSource, // [kg.m/s] (tang'l momentum)
const scalar pressureSource, // [kg.m/s] (normal momentum)
const scalar energySource = 0 // [J]
);
//- Add sources
virtual void addSources
(
const label patchi, // patchi on primary region
const label facei, // facei of patchi
const scalar massSource, // [kg]
const vector& momentumSource, // [kg.m/s] (tang'l momentum)
const scalar pressureSource, // [kg.m/s] (normal momentum)
const scalar energySource = 0 // [J]
);
};

View File

@ -31,7 +31,6 @@ License
#include "gravityMeshObject.H"
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
@ -91,9 +90,7 @@ liquidFilmModel::liquidFilmModel
)
:
liquidFilmBase(modelType, mesh, dict),
thermo_(dict.subDict("thermo")),
rho_
(
IOobject
@ -221,7 +218,6 @@ liquidFilmModel::liquidFilmModel
dimensionedScalar(dimMass, Zero),
calculatedFvPatchField<scalar>::typeName
),
cloudDiameterTrans_
(
IOobject
@ -236,79 +232,22 @@ liquidFilmModel::liquidFilmModel
dimensionedScalar(dimLength, Zero),
calculatedFvPatchField<scalar>::typeName
),
turbulence_(filmTurbulenceModel::New(*this, dict)),
turbulencePtr_(filmTurbulenceModel::New(*this, dict)),
availableMass_(regionMesh().faces().size(), Zero),
injection_(*this, dict),
forces_(*this, dict)
{
if (dict.readIfPresent("T0", Tref_))
{
Tf_ = dimensionedScalar("T0", dimTemperature, dict);
}
correctThermoFields();
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
const areaScalarField& liquidFilmModel::mu() const
{
return mu_;
}
const areaScalarField& liquidFilmModel::rho() const
{
return rho_;
}
const areaScalarField& liquidFilmModel::sigma() const
{
return sigma_;
}
const areaScalarField& liquidFilmModel::Tf() const
{
return Tf_;
}
const areaScalarField& liquidFilmModel::Cp() const
{
return Cp_;
}
const liquidMixtureProperties& liquidFilmModel::thermo() const
{
return thermo_;
}
scalar liquidFilmModel::Tref() const
{
return Tref_;
}
const volScalarField& liquidFilmModel::cloudMassTrans() const
{
return cloudMassTrans_;
}
const volScalarField& liquidFilmModel::cloudDiameterTrans() const
{
return cloudDiameterTrans_;
}
void liquidFilmModel::preEvolveRegion()
{
liquidFilmBase::preEvolveRegion();
@ -361,12 +300,11 @@ void liquidFilmModel::info()
const DimensionedField<scalar, areaMesh>& sf = regionMesh().S();
Info<< indent << "min/max(mag(Uf)) = "
<< gMin(mag(Uf_.field())) << ", "
<< gMax(mag(Uf_.field())) << nl
<< gMinMaxMag(Uf_.field()) << nl
<< indent << "min/max(delta) = "
<< gMin(h_.field()) << ", " << gMax(h_.field()) << nl
<< gMinMax(h_.field()) << nl
<< indent << "coverage = "
<< gSum(alpha()().field()*mag(sf.field()))/gSum(mag(sf.field())) << nl
<< gSum(alpha()().field()*mag(sf.field()))/gSumMag(sf.field()) << nl
<< indent << "total mass = "
<< gSum(availableMass_) << nl;

View File

@ -27,17 +27,44 @@ Class
Foam::regionFaModels::liquidFilmModel
Description
Thin Model Film model.
Thin film model.
Usage
Example of the boundary condition specification:
\verbatim
<patchName>
{
// Mandatory entries
thermo <dict>;
// Optional entries
T0 <scalar>;
// Inherited entries
...
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
thermo | Liquid thermodynamics properties | dict | yes | -
T0 | Reference temperature [K] | scalar | choice | -
\endtable
The inherited entries are elaborated in:
- \link liquidFilmBase.H \endlink
- \link filmTurbulenceModel.H \endlink
- \link injectionModelList.H \endlink
- \link forceList.H \endlink
SourceFiles
liquidFilmModel.C
kinematicThinFilmI.H
\*---------------------------------------------------------------------------*/
#ifndef liquidFilmModel_H
#define liquidFilmModel_H
#ifndef areaSurfaceFilmModels_liquidFilmModel_H
#define areaSurfaceFilmModels_liquidFilmModel_H
#include "volFieldsFwd.H"
#include "liquidFilmBase.H"
@ -58,7 +85,7 @@ namespace areaSurfaceFilmModels
{
/*---------------------------------------------------------------------------*\
Class liquidFilmModel Declaration
Class liquidFilmModel Declaration
\*---------------------------------------------------------------------------*/
class liquidFilmModel
@ -72,7 +99,7 @@ protected:
//- Liquid thermo
liquidMixtureProperties thermo_;
//- Reference tempararure
//- Reference temperature
scalar Tref_;
@ -84,16 +111,16 @@ protected:
//- Dynamic viscosity [Pa.s]
areaScalarField mu_;
//- Film temperature
//- Film temperature [K]
areaScalarField Tf_;
//- Film Heat capacity
//- Film heat capacity [J/K]
areaScalarField Cp_;
//- Surface tension [m/s2]
//- Surface tension [kg/s^2]
areaScalarField sigma_;
//- Film rho*height
//- Film rho*height [m.kg/m^3]
areaScalarField hRho_;
@ -118,25 +145,19 @@ protected:
volScalarField cloudDiameterTrans_;
// General properties
// Submodels
//- Turbulence model
autoPtr<filmTurbulenceModel> turbulence_;
autoPtr<filmTurbulenceModel> turbulencePtr_;
//- Available mass for transfer via sub-models
scalarField availableMass_;
// Sub-models
//- Cloud injection
injectionModelList injection_;
//- Available mass for transfer via sub-models
scalarField availableMass_;
//- Cloud injection
injectionModelList injection_;
//- Transfer with the continuous phase
//transferModelList transfer_;
//- List of film forces
forceList forces_;
//- List of film forces
forceList forces_;
public:
@ -168,58 +189,67 @@ public:
// Member Functions
// Helpers
// Helpers
//- Correct thermo
void correctThermoFields();
//- Correct thermo
void correctThermoFields();
// Access
// Access
//- Access const reference mu
const areaScalarField& mu() const;
//- Return const reference mu
const areaScalarField& mu() const noexcept { return mu_; }
//- Access const reference rho
const areaScalarField& rho() const;
//- Return const reference rho
const areaScalarField& rho() const noexcept { return rho_; }
//- Access const reference sigma
const areaScalarField& sigma() const;
//- Return const reference sigma
const areaScalarField& sigma() const noexcept { return sigma_; }
//- Access const reference Tf
const areaScalarField& Tf() const;
//- Return const reference Tf
const areaScalarField& Tf() const noexcept { return Tf_; }
//- Access const reference Cp
const areaScalarField& Cp() const;
//- Return const reference Cp
const areaScalarField& Cp() const noexcept { return Cp_; }
//- Access to thermo
const liquidMixtureProperties& thermo() const;
//- Return const reference thermo
const liquidMixtureProperties& thermo() const noexcept
{
return thermo_;
}
//- Access to reference temperature
scalar Tref() const;
//- Return reference temperature
scalar Tref() const noexcept { return Tref_; }
// Transfer fields - to the primary region (lagragian injection)
// Transfer fields - to the primary region (lagragian injection)
//- Return the film mass available for transfer to cloud
virtual const volScalarField& cloudMassTrans() const;
//- Return the film mass available for transfer to cloud
virtual const volScalarField& cloudMassTrans() const noexcept
{
return cloudMassTrans_;
}
//- Return the parcel diameters originating from film to cloud
virtual const volScalarField& cloudDiameterTrans() const;
//- Return the parcel diameters originating from film to cloud
virtual const volScalarField& cloudDiameterTrans() const noexcept
{
return cloudDiameterTrans_;
}
// Evolution
// Evolution
//- Pre-evolve film
virtual void preEvolveRegion();
//- Pre-evolve film
virtual void preEvolveRegion();
//- Post-evolve film
virtual void postEvolveRegion();
//- Post-evolve film
virtual void postEvolveRegion();
// I-O
// I-O
//- Provide some feedback
virtual void info();
//- Provide some feedback
virtual void info();
};

View File

@ -35,8 +35,8 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef filmSubModelBase_H
#define filmSubModelBase_H
#ifndef areaSurfaceFilmModels_filmSubModelBase_H
#define areaSurfaceFilmModels_filmSubModelBase_H
#include "liquidFilmBase.H"
#include "subModelBase.H"
@ -100,19 +100,19 @@ public:
// Member Functions
// Access
// Access
//- Flag to indicate when to write a property
virtual bool writeTime() const;
//- Flag to indicate when to write a property
virtual bool writeTime() const;
//- Return const access to the film surface film model
inline const liquidFilmBase& film() const;
//- Return const access to the film surface film model
inline const liquidFilmBase& film() const;
//- Return the reference to the film surface film model
inline liquidFilmBase& film();
//- Return the reference to the film surface film model
inline liquidFilmBase& film();
template<class FilmType>
inline const FilmType& filmType() const;
template<class FilmType>
inline const FilmType& filmType() const;
};

View File

@ -94,26 +94,19 @@ filmTurbulenceModel::filmTurbulenceModel
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
const liquidFilmBase& filmTurbulenceModel::film() const
{
return film_;
}
tmp<areaScalarField> filmTurbulenceModel::Cw() const
{
auto tCw =
tmp<areaScalarField>::New
auto tCw = tmp<areaScalarField>::New
(
IOobject
(
IOobject
(
"tCw",
film_.primaryMesh().time().timeName(),
film_.primaryMesh()
),
film_.regionMesh(),
dimensionedScalar(dimVelocity)
);
"tCw",
film_.primaryMesh().time().timeName(),
film_.primaryMesh()
),
film_.regionMesh(),
dimensionedScalar(dimVelocity)
);
auto& Cw = tCw.ref();
@ -195,9 +188,9 @@ tmp<faVectorMatrix> filmTurbulenceModel::primaryRegionFriction
areaVectorField& U
) const
{
tmp<faVectorMatrix> tshearStress
auto tshearStress = tmp<faVectorMatrix>::New
(
new faVectorMatrix(U, sqr(U.dimensions())*sqr(dimLength))
U, sqr(U.dimensions())*sqr(dimLength)
);
switch (shearMethod_)

View File

@ -29,14 +29,77 @@ Class
Description
Base class for film turbulence models
Usage
Example of the model specification:
\verbatim
{
// Mandatory entries
turbulence <model>;
<model>Coeffs
{
// Mandatory entries
friction <word>;
shearStress <word>;
// Optional entries
rho <word>;
// Conditional entries
// if rho=rhoInf
rhoInf <scalar>;
// if friction=DarcyWeisbach
DarcyWeisbach <scalar>;
// if friction=ManningStrickler
n <scalar>;
// if shearStress=simple
Cf <scalar>;
}
// Inherited entries
...
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
turbulence | Model name | word | yes | -
friction | Friction model | word | yes | -
shearStress | Shear-stress model | word | yes | -
rho | Name of density field | word | no | rho
rhoInf | Reference density | scalar | choice | -
DarcyWeisbach | Friction-model constant | scalar | choice | -
n | Friction-model constant | scalar | choice | -
Cf | Skin-friction coefficient | scalar | choice | -
\endtable
Options for the \c friction entry:
\verbatim
quadraticProfile | Quadratic-profile model
linearProfile | Linear-profile model
DarcyWeisbach | Darcy-Weisbach model
ManningStrickler | Manning-Strickler model
\endverbatim
Options for the \c shearStress entry:
\verbatim
simple | Skin-friction coefficient model
wallFunction | Wall-function model
\endverbatim
SourceFiles
filmTurbulenceModel.C
filmTurbulenceModelNew.C
\*---------------------------------------------------------------------------*/
#ifndef regionFaModels_filmTurbulenceModel_H
#define regionFaModels_filmTurbulenceModel_H
#ifndef areaSurfaceFilmModels_filmTurbulenceModel_H
#define areaSurfaceFilmModels_filmTurbulenceModel_H
#include "areaFieldsFwd.H"
#include "runTimeSelectionTables.H"
@ -106,13 +169,13 @@ protected:
//- Model dictionary
const dictionary dict_;
//- Method used
//- Friction model
const frictionMethodType method_;
//- Shear method used
//- Shear-stress model
const shearMethodType shearMethod_;
//- Name of density field (optional)
//- Name of density field
word rhoName_;
//- Reference density needed for incompressible calculations
@ -167,40 +230,37 @@ public:
// Member Functions
// Access
// Access
//- Return film
const liquidFilmBase& film() const;
//- Return const reference to film
const liquidFilmBase& film() const { return film_; }
// Turbulence
// Turbulence
//- Return the effective viscous stress (laminar + turbulent)
tmp<volSymmTensorField> devRhoReff() const;
//- Return the effective viscous stress (laminar + turbulent)
tmp<volSymmTensorField> devRhoReff() const;
//- Return primary region friction
tmp<faVectorMatrix> primaryRegionFriction
(
areaVectorField& U
) const;
//- Return primary region friction
tmp<faVectorMatrix> primaryRegionFriction
(
areaVectorField& U
) const;
//- Return rho if specified otherwise rhoRef
tmp<volScalarField> rho() const;
//- Return rho if specified otherwise rhoRef
tmp<volScalarField> rho() const;
//- Return the wall film surface friction
virtual tmp<areaScalarField> Cw() const;
//- Return the film turbulence viscosity
virtual tmp<areaScalarField> mut() const = 0;
//- Return the wall film surface friction
virtual tmp<areaScalarField> Cw() const;
// Evaluation
// Evaluation
//- Correct/update the model
virtual void correct() = 0;
//- Correct/update the model
virtual void correct() = 0;
//- Return the source for the film momentum equation
virtual tmp<faVectorMatrix> Su(areaVectorField& U) const = 0;
//- Return the source for the film momentum equation
virtual tmp<faVectorMatrix> Su(areaVectorField& U) const = 0;
};

View File

@ -28,7 +28,6 @@ License
#include "laminar.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
@ -43,6 +42,22 @@ namespace areaSurfaceFilmModels
defineTypeNameAndDebug(laminar, 0);
addToRunTimeSelectionTable(filmTurbulenceModel, laminar, dictionary);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
tmp<faVectorMatrix> laminar::wallFriction(areaVectorField& U) const
{
// Local references to film fields
tmp<areaVectorField> Uw = film_.Uw();
tmp<areaScalarField> wf = Cw();
return
(
- fam::Sp(wf(), U) + wf()*Uw() // wall contribution
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
laminar::laminar
@ -57,26 +72,6 @@ laminar::laminar
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
tmp<areaScalarField> laminar::mut() const
{
auto tmut = tmp<areaScalarField>::New
(
IOobject
(
"mut",
film().primaryMesh().time().timeName(),
film().primaryMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
film().regionMesh(),
dimensionedScalar(dimMass/dimLength/dimTime)
);
return tmut;
}
void laminar::correct()
{}
@ -87,19 +82,6 @@ tmp<faVectorMatrix> laminar::Su(areaVectorField& U) const
}
tmp<faVectorMatrix> laminar::wallFriction(areaVectorField& U) const
{
// local references to film fields
tmp<areaVectorField> Uw = film_.Uw();
tmp<areaScalarField> wf = Cw();
return
(
- fam::Sp(wf(), U) + wf()*Uw() // wall contribution
);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace areaSurfaceFilmModels

View File

@ -34,10 +34,10 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef regionFaModels_surfaceFilmModels_laminar_H
#define regionFaModels_surfaceFilmModels_laminar_H
#ifndef areaSurfaceFilmModels_surfaceFilmModels_laminar_H
#define areaSurfaceFilmModels_surfaceFilmModels_laminar_H
#include "liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H"
#include "filmTurbulenceModel.H"
#include "faCFD.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -59,6 +59,9 @@ class laminar
{
// Private Member Functions
//- Return wall-friction system
tmp<faVectorMatrix> wallFriction(areaVectorField& U) const;
//- No copy construct
laminar(const laminar&) = delete;
@ -84,22 +87,16 @@ public:
// Member Functions
// Turbulence
// Turbulence
//- Wall friction
tmp<faVectorMatrix> wallFriction(areaVectorField& U) const;
//- Return the film turbulence viscosity
virtual tmp<areaScalarField> mut() const;
//- Return the source for the film momentum equation
virtual tmp<faVectorMatrix> Su(areaVectorField& U) const;
//- Return the source for the film momentum equation
virtual tmp<faVectorMatrix> Su(areaVectorField& U) const;
// Evaluation
// Evaluation
//- Correct/update the model
virtual void correct();
//- Correct/update the model
virtual void correct();
};

View File

@ -35,8 +35,8 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef force_H
#define force_H
#ifndef areaSurfaceFilmModels_force_H
#define areaSurfaceFilmModels_force_H
#include "filmSubModelBase.H"
#include "runTimeSelectionTables.H"
@ -87,6 +87,7 @@ public:
(film, dict)
);
// Constructors
//- Construct null
@ -118,10 +119,10 @@ public:
// Member Functions
// Evolution
// Evaluation
//- Correct
virtual tmp<faVectorMatrix> correct(areaVectorField& U) = 0;
//- Correct
virtual tmp<faVectorMatrix> correct(areaVectorField& U) = 0;
};

Some files were not shown because too many files have changed in this diff Show More