zhyang-dev/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ENH: add coupledFaPatch::delta()

Browse Source 
- symmetrical evaluation for processor patches, eliminates
  scalar/vector multiply followed by projection.

STYLE: use evaluateCoupled instead of local versions
This commit is contained in:
Mark Olesen
2023-01-26 17:12:45 +01:00
parent eb8b51b475
commit ea2bedf073
9 changed files with 80 additions and 121 deletions
Download Patch File Download Diff File Expand all files Collapse all files

49
applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C
View File

@ -50,8 +50,7 @@ Description
#include "turbulentTransportModel.H"
#include "turbulentFluidThermoModel.H"
#include "wallDist.H"
#include "processorFvPatchField.H"
#include "zeroGradientFvPatchField.H"
#include "processorFvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -61,35 +60,11 @@ static const scalar kappa(0.41);
template<class Type>
void correctProcessorPatches
(
GeometricField<Type, fvPatchField, volMesh>& vf
)
void correctProcessorPatches(GeometricField<Type, fvPatchField, volMesh>& fld)
{
if (!Pstream::parRun())
if (UPstream::parRun())
{
return;
}
// Not possible to use correctBoundaryConditions on fields as they may
// use local info as opposed to the constraint values employed here,
// but still need to update processor patches
auto& bf = vf.boundaryFieldRef();
forAll(bf, patchi)
{
if (isA<processorFvPatchField<Type>>(bf[patchi]))
{
bf[patchi].initEvaluate();
}
}
forAll(bf, patchi)
{
if (isA<processorFvPatchField<Type>>(bf[patchi]))
{
bf[patchi].evaluate();
}
fld.boundaryFieldRef().template evaluateCoupled<processorFvPatch>();
}
}
@ -119,11 +94,11 @@ void blendField
pf = (1 - mask)*pf + mask*boundaryLayerField;
fld.max(SMALL);
// Correct the processor patches only.
// Do not correct BC
// - operation may use inconsistent fields wrt these local
// manipulations
//fld.correctBoundaryConditions();
correctProcessorPatches<scalar>(fld);
correctProcessorPatches(fld);
Info<< "Writing " << fieldName << nl << endl;
fld.write();
@ -158,11 +133,11 @@ void calcOmegaField
pf = (1 - mask)*pf + mask*epsilonBL/(Cmu*kBL + SMALL);
omega.max(SMALL);
// Correct the processor patches only.
// Do not correct BC
// - operation may use inconsistent fields wrt these local
// manipulations
// omega.correctBoundaryConditions();
correctProcessorPatches<scalar>(omega);
correctProcessorPatches(omega);
Info<< "Writing omega\n" << endl;
omega.write();
@ -192,11 +167,11 @@ void setField
volScalarField fld(fldHeader, mesh);
fld = value;
// Correct the processor patches only.
// Do not correct BC
// - operation may use inconsistent fields wrt these local
// manipulations
// fld.correctBoundaryConditions();
correctProcessorPatches<scalar>(fld);
correctProcessorPatches(fld);
Info<< "Writing " << fieldName << nl << endl;
fld.write();
@ -343,7 +318,7 @@ int main(int argc, char *argv[])
}
}
mask.correctBoundaryConditions();
correctProcessorPatches<vector>(U);
correctProcessorPatches(U);
if (writeTurbulenceFields)
{
@ -356,7 +331,7 @@ int main(int argc, char *argv[])
// Do not correct BC - wall functions will 'undo' manipulation above
// by using nut from turbulence model
correctProcessorPatches<scalar>(nut);
correctProcessorPatches(nut);
Info<< "Writing nut\n" << endl;
nut.write();

46
applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C
View File

@ -109,8 +109,8 @@ Note
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "turbulentFluidThermoModel.H"
#include "processorFvPatchField.H"
#include "wallFvPatch.H"
#include "processorFvPatch.H"
#include "fixedValueFvPatchFields.H"
using namespace Foam;
@ -124,35 +124,11 @@ void InfoField(const word& fldName)
template<class Type>
void correctProcessorPatches
(
GeometricField<Type, fvPatchField, volMesh>& vf
)
void correctProcessorPatches(GeometricField<Type, fvPatchField, volMesh>& fld)
{
if (!Pstream::parRun())
if (UPstream::parRun())
{
return;
}
// Not possible to use correctBoundaryConditions on fields as they may
// use local info as opposed to the constraint values employed here,
// but still need to update processor patches
auto& bf = vf.boundaryFieldRef();
forAll(bf, patchi)
{
if (isA<processorFvPatchField<Type>>(bf[patchi]))
{
bf[patchi].initEvaluate();
}
}
forAll(bf, patchi)
{
if (isA<processorFvPatchField<Type>>(bf[patchi]))
{
bf[patchi].evaluate();
}
fld.boundaryFieldRef().template evaluateCoupled<processorFvPatch>();
}
}
@ -411,7 +387,7 @@ int main(int argc, char *argv[])
),
mesh,
dimensionedScalar(dimless, scalar(1)),
fixedValueFvPatchField<scalar>::typeName
fixedValueFvPatchScalarField::typeName
);
for (fvPatchScalarField& pfld : f.boundaryFieldRef())
@ -482,31 +458,31 @@ int main(int argc, char *argv[])
// (M:Eq. 9)
const dimensionedScalar maxU(dimVelocity, SMALL);
U *= min(scalar(1), fRei*uTau/max(mag(U), maxU));
correctProcessorPatches<vector>(U);
correctProcessorPatches(U);
}
if (tepsilon.valid())
{
tepsilon.ref() = epsilon;
correctProcessorPatches<scalar>(tepsilon.ref());
correctProcessorPatches(tepsilon.ref());
}
if (tk.valid())
{
tk.ref() = k;
correctProcessorPatches<scalar>(tk.ref());
correctProcessorPatches(tk.ref());
}
if (tomega.valid())
{
const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
tomega.ref() = Cmu*epsilon/(k + k0);
correctProcessorPatches<scalar>(tomega.ref());
correctProcessorPatches(tomega.ref());
}
if (tR.valid())
{
volSymmTensorField& R = tR.ref();
auto& R = tR.ref();
// (M:Eq. 3)
const volSphericalTensorField Rdiag(k*twoThirdsI);
@ -514,7 +490,7 @@ int main(int argc, char *argv[])
{
R[celli] = Rdiag[celli];
}
correctProcessorPatches<symmTensor>(R);
correctProcessorPatches(R);
}

46
src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.C
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -32,7 +33,7 @@ License
namespace Foam
{
defineTypeNameAndDebug(coupledFaPatch, 0);
defineTypeName(coupledFaPatch);
}
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
@ -47,25 +48,28 @@ void Foam::coupledFaPatch::calcTransformTensors
{
if (mag(nf & nr) < 1 - SMALL)
{
separation_.setSize(0);
separation_.clear();
forwardT_.resize(1);
reverseT_.resize(1);
forwardT_ = tensorField(1, rotationTensor(-nr, nf));
reverseT_ = tensorField(1, rotationTensor(nf, -nr));
forwardT_[0] = rotationTensor(-nr, nf);
reverseT_[0] = rotationTensor(nf, -nr);
}
else
{
forwardT_.setSize(0);
reverseT_.setSize(0);
forwardT_.clear();
reverseT_.clear();
vector separation = (nf & (Cr - Cf))*nf;
if (mag(separation) > SMALL)
{
separation_ = vectorField(1, separation);
separation_.resize(1);
separation_[0] = separation;
}
else
{
separation_.setSize(0);
separation_.clear();
}
}
}
@ -81,10 +85,10 @@ void Foam::coupledFaPatch::calcTransformTensors
{
if (sum(mag(nf & nr)) < Cf.size() - SMALL)
{
separation_.setSize(0);
separation_.clear();
forwardT_.setSize(size());
reverseT_.setSize(size());
forwardT_.resize_nocopy(size());
reverseT_.resize_nocopy(size());
forAll(forwardT_, facei)
{
@ -94,27 +98,35 @@ void Foam::coupledFaPatch::calcTransformTensors
if (sum(mag(forwardT_ - forwardT_[0])) < SMALL)
{
forwardT_.setSize(1);
reverseT_.setSize(1);
forwardT_.resize(1);
reverseT_.resize(1);
}
}
else
{
forwardT_.setSize(0);
reverseT_.setSize(0);
forwardT_.clear();
reverseT_.clear();
separation_ = (nf&(Cr - Cf))*nf;
if (sum(mag(separation_)) < SMALL)
{
separation_.setSize(0);
separation_.clear();
}
else if (sum(mag(separation_ - separation_[0])) < SMALL)
{
separation_.setSize(1);
separation_.resize(1);
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::vectorField> Foam::coupledFaPatch::delta() const
{
return (edgeCentres() - edgeFaceCentres());
}
// ************************************************************************* //

2
src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H
View File

@ -104,7 +104,7 @@ protected:
public:
//- Runtime type information
TypeName("coupled");
TypeNameNoDebug("coupled");
// Constructors

10
src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C
View File

@ -169,7 +169,7 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const
{
const scalarField& magL = magEdgeLengths();
const scalarField deltas(edgeNormals() & faPatch::delta());
const scalarField deltas(edgeNormals() & coupledFaPatch::delta());
const label sizeby2 = deltas.size()/2;
scalar maxMatchError = 0;
@ -222,7 +222,7 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const
void Foam::cyclicFaPatch::makeDeltaCoeffs(scalarField& dc) const
{
const scalarField deltas(edgeNormals() & faPatch::delta());
const scalarField deltas(edgeNormals() & coupledFaPatch::delta());
const label sizeby2 = deltas.size()/2;
for (label edgei = 0; edgei < sizeby2; ++edgei)
@ -272,7 +272,7 @@ void Foam::cyclicFaPatch::movePoints
Foam::tmp<Foam::vectorField> Foam::cyclicFaPatch::delta() const
{
const vectorField patchD(faPatch::delta());
const vectorField patchD(coupledFaPatch::delta());
const label sizeby2 = patchD.size()/2;
auto tpdv = tmp<vectorField>::New(patchD.size());
@ -284,7 +284,7 @@ Foam::tmp<Foam::vectorField> Foam::cyclicFaPatch::delta() const
for (label edgei = 0; edgei < sizeby2; ++edgei)
{
const vector& ddi = patchD[edgei];
vector dni = patchD[edgei + sizeby2];
const vector& dni = patchD[edgei + sizeby2];
pdv[edgei] = ddi - dni;
pdv[edgei + sizeby2] = -pdv[edgei];
@ -295,7 +295,7 @@ Foam::tmp<Foam::vectorField> Foam::cyclicFaPatch::delta() const
for (label edgei = 0; edgei < sizeby2; ++edgei)
{
const vector& ddi = patchD[edgei];
vector dni = patchD[edgei + sizeby2];
const vector& dni = patchD[edgei + sizeby2];
pdv[edgei] = ddi - transform(forwardT()[0], dni);
pdv[edgei + sizeby2] = -transform(reverseT()[0], pdv[edgei]);

24
src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2019-2022 OpenCFD Ltd.
Copyright (C) 2019-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -395,7 +395,7 @@ void Foam::processorFaPatch::makeWeights(scalarField& w) const
(
(
neighbEdgeLengths()
/mag(neighbEdgeLengths())
/(mag(neighbEdgeLengths()) + VSMALL)
)
& (
neighbEdgeCentres()
@ -405,13 +405,13 @@ void Foam::processorFaPatch::makeWeights(scalarField& w) const
w = neighbEdgeCentresCn/
(
(edgeNormals() & faPatch::delta())
+ neighbEdgeCentresCn
(edgeNormals() & coupledFaPatch::delta())
+ neighbEdgeCentresCn + VSMALL
);
}
else
{
w = 1.0;
w = scalar(1);
}
}
@ -420,11 +420,13 @@ void Foam::processorFaPatch::makeDeltaCoeffs(scalarField& dc) const
{
if (Pstream::parRun())
{
dc = (1.0 - weights())/(edgeNormals() & faPatch::delta());
dc = (1.0 - weights())
/((edgeNormals() & coupledFaPatch::delta()) + VSMALL);
}
else
{
dc = 1.0/(edgeNormals() & faPatch::delta());
dc = scalar(1)
/((edgeNormals() & coupledFaPatch::delta()) + VSMALL);
}
}
@ -433,11 +435,11 @@ Foam::tmp<Foam::vectorField> Foam::processorFaPatch::delta() const
{
if (Pstream::parRun())
{
// To the transformation if necessary
// Do the transformation if necessary
if (parallel())
{
return
faPatch::delta()
coupledFaPatch::delta()
- (
neighbEdgeCentres()
- neighbEdgeFaceCentres()
@ -446,7 +448,7 @@ Foam::tmp<Foam::vectorField> Foam::processorFaPatch::delta() const
else
{
return
faPatch::delta()
coupledFaPatch::delta()
- transform
(
forwardT(),
@ -459,7 +461,7 @@ Foam::tmp<Foam::vectorField> Foam::processorFaPatch::delta() const
}
else
{
return faPatch::delta();
return coupledFaPatch::delta();
}
}

6
src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
View File

@ -482,11 +482,11 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::delta() const
vectorField edgePN(edgeCentres() - edgeFaceCentres());
// Do not allow any mag(val) < SMALL
for (vector& edgei : edgePN)
for (vector& e : edgePN)
{
if (mag(edgei) < SMALL)
if (e.magSqr() < ROOTSMALL)
{
edgei = vector::uniform(SMALL);
e = vector::uniform(SMALL);
}
}

8
src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.C
View File

@ -31,16 +31,10 @@ License
namespace Foam
{
defineTypeNameAndDebug(coupledFvPatch, 0);
defineTypeName(coupledFvPatch);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::coupledFvPatch::~coupledFvPatch()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::vectorField> Foam::coupledFvPatch::delta() const

10
src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H
View File

@ -35,8 +35,8 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef coupledFvPatch_H
#define coupledFvPatch_H
#ifndef Foam_coupledFvPatch_H
#define Foam_coupledFvPatch_H
#include "fvPatch.H"
#include "lduInterface.H"
@ -56,7 +56,7 @@ class coupledFvPatch
public lduInterface,
public fvPatch
{
// Private data
// Private Data
const coupledPolyPatch& coupledPolyPatch_;
@ -75,7 +75,7 @@ public:
//- Runtime type information
TypeName(coupledPolyPatch::typeName_());
TypeNameNoDebug(coupledPolyPatch::typeName_());
// Constructors
@ -89,7 +89,7 @@ public:
//- Destructor
virtual ~coupledFvPatch();
virtual ~coupledFvPatch() = default;
// Member Functions