From 2d731e1af338d038ebbb310c19da7f6394f7f91b Mon Sep 17 00:00:00 2001 From: Kutalmis Bercin Date: Fri, 1 Nov 2024 13:56:35 +0000 Subject: [PATCH] ENH: processorFaPatch: use internal-edge algos for processor edges to ensure parallel consistency - The edgeInterpolation::makeCorrectionVectors() disables the non-orthogonality correction if the calculated non-orthogonality coefficient is below 0.1. However, this activation routine only considers internal edges, and excludes any processor edges, resulting in inconsistent parallel calculations. This routine is removed. - Fatal errors are replaced with zero-valued fields for non-orthogonality- and skewness-correction routines. --- .../faPatches/basic/coupled/coupledFaPatch.H | 6 + .../constraint/cyclic/cyclicFaPatch.C | 7 ++ .../constraint/cyclic/cyclicFaPatch.H | 6 + .../constraint/processor/processorFaPatch.C | 112 ++++++++++++++++-- .../constraint/processor/processorFaPatch.H | 6 + .../faMesh/faPatches/faPatch/faPatch.C | 23 +++- .../faMesh/faPatches/faPatch/faPatch.H | 14 ++- .../edgeInterpolation/edgeInterpolation.C | 58 ++++----- .../system/finite-area/faSchemes | 2 +- 9 files changed, 184 insertions(+), 50 deletions(-) diff --git a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H index a4bdba28e5..d7e0dc74ce 100644 --- a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H +++ b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H @@ -79,9 +79,15 @@ protected: //- Make patch weighting factors virtual void makeWeights(scalarField&) const = 0; + //- Make patch geodesic distance between P and N + virtual void makeLPN(scalarField&) const = 0; + //- Make patch face - neighbour cell distances virtual void makeDeltaCoeffs(scalarField&) const = 0; + //- Make non-orthogonality correction vectors + virtual void makeCorrectionVectors(vectorField&) const = 0; + //- Calculate the uniform transformation tensors void calcTransformTensors ( diff --git a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C index 5d0ffcb436..426cc28158 100644 --- a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C +++ b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C @@ -218,6 +218,13 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const } +void Foam::cyclicFaPatch::makeLPN(scalarField& lPN) const +{ + makeDeltaCoeffs(lPN); + lPN = scalar(1)/lPN; +} + + void Foam::cyclicFaPatch::makeDeltaCoeffs(scalarField& dc) const { const scalarField deltas(edgeNormals() & coupledFaPatch::delta()); diff --git a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H index 9e4059fdff..e9c559353a 100644 --- a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H +++ b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H @@ -78,9 +78,15 @@ protected: //- Make patch weighting factors void makeWeights(scalarField&) const; + //- Make patch geodesic distance between P and N + void makeLPN(scalarField&) const; + //- Make patch face - neighbour cell distances void makeDeltaCoeffs(scalarField&) const; + //- Make non-orthogonality correction vectors + void makeCorrectionVectors(vectorField& cv) const { cv = Zero; } + public: diff --git a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C index c60a08a0ba..d38a88c165 100644 --- a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C +++ b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C @@ -398,18 +398,21 @@ void Foam::processorFaPatch::makeWeights(scalarField& w) const { if (Pstream::parRun()) { - // The face normals point in the opposite direction on the other side - scalarField neighbEdgeCentresCn + const vectorField& skew = skewCorrectionVectors(); + const scalarField& PN = lPN(); + + const scalarField lEN ( - neighbEdgeNormals() - & (neighbEdgeCentres() - neighbEdgeFaceCentres()) + mag(neighbEdgeFaceCentres() - edgeCentres() + skew) ); - w = neighbEdgeCentresCn/ - ( - (edgeNormals() & coupledFaPatch::delta()) - + neighbEdgeCentresCn + VSMALL - ); + forAll(w, i) + { + if (mag(PN[i]) > SMALL) + { + w[i] = lEN[i]/PN[i]; + } + } } else { @@ -418,12 +421,61 @@ void Foam::processorFaPatch::makeWeights(scalarField& w) const } +void Foam::processorFaPatch::makeLPN(scalarField& lPN) const +{ + const vectorField& skew = skewCorrectionVectors(); + + if (Pstream::parRun()) + { + lPN = + mag(edgeCentres() - skew - edgeFaceCentres()) // lPE + + mag(neighbEdgeFaceCentres() - edgeCentres() + skew); // lEN + + for (scalar& lpn: lPN) + { + if (mag(lpn) < SMALL) + { + lpn = SMALL; + } + } + } + else + { + lPN = + mag(edgeCentres() - skew - edgeFaceCentres()) // lPE + + mag(edgeFaceCentres() - edgeCentres() + skew); // lEN + } +} + + void Foam::processorFaPatch::makeDeltaCoeffs(scalarField& dc) const { if (Pstream::parRun()) { - dc = (1.0 - weights()) - /((edgeNormals() & coupledFaPatch::delta()) + VSMALL); + const edgeList& edges = boundaryMesh().mesh().edges(); + const pointField& points = boundaryMesh().mesh().points(); + const vectorField& lengths = edgeLengths(); + const scalarField& PN = lPN(); + + tmp tdelta = processorFaPatch::delta(); + vectorField& unitDelta = tdelta.ref(); + + forAll(dc, i) + { + vector edgeNormal = + normalised(lengths[i] ^ edges[i + start()].vec(points)); + + unitDelta[i].removeCollinear(edgeNormal); + unitDelta[i].normalise(); + + edgeNormal = normalised(lengths[i]); + + const scalar alpha = PN[i]*(edgeNormal & unitDelta[i]); + if (mag(alpha) > SMALL) + { + dc[i] = scalar(1)/max(alpha, 0.05*PN[i]); + } + } } else { @@ -433,6 +485,44 @@ void Foam::processorFaPatch::makeDeltaCoeffs(scalarField& dc) const } +void Foam::processorFaPatch::makeCorrectionVectors(vectorField& cv) const +{ + if (Pstream::parRun()) + { + const edgeList& edges = boundaryMesh().mesh().edges(); + const pointField& points = boundaryMesh().mesh().points(); + const vectorField& lengths = edgeLengths(); + + tmp tdelta = processorFaPatch::delta(); + vectorField& unitDelta = tdelta.ref(); + + forAll(cv, i) + { + vector edgeNormal = + normalised(lengths[i] ^ edges[i + start()].vec(points)); + + unitDelta[i].removeCollinear(edgeNormal); + unitDelta.normalise(); + + edgeNormal = normalised(lengths[i]); + + const scalar alpha = unitDelta[i] & edgeNormal; + scalar dc = SMALL; + if (mag(alpha) > SMALL) + { + dc = scalar(1)/alpha; + } + + cv[i] = edgeNormal - dc*unitDelta[i]; + } + } + else + { + cv = Zero; + } +} + + Foam::tmp Foam::processorFaPatch::delta() const { if (Pstream::parRun()) diff --git a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H index 6bcf0b05f9..7eeb4eb3aa 100644 --- a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H +++ b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H @@ -105,9 +105,15 @@ protected: //- Make patch weighting factors void makeWeights(scalarField&) const; + //- Make patch geodesic distance between P and N + void makeLPN(scalarField&) const; + //- Make patch face - neighbour cell distances void makeDeltaCoeffs(scalarField&) const; + //- Make non-orthogonality correction vectors + void makeCorrectionVectors(vectorField&) const; + //- Find non-globa patch points void makeNonGlobalPatchPoints() const; diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C index 735b9b2224..9b103bd872 100644 --- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C +++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C @@ -500,12 +500,30 @@ Foam::tmp Foam::faPatch::delta() const } +void Foam::faPatch::makeLPN(scalarField& lPN) const +{ + lPN = (edgeNormals() & delta()); +} + + +const Foam::scalarField& Foam::faPatch::lPN() const +{ + return boundaryMesh().mesh().lPN().boundaryField()[index()]; +} + + void Foam::faPatch::makeDeltaCoeffs(scalarField& dc) const { dc = scalar(1)/(edgeNormals() & delta()); } +const Foam::scalarField& Foam::faPatch::deltaCoeffs() const +{ + return boundaryMesh().mesh().deltaCoeffs().boundaryField()[index()]; +} + + void Foam::faPatch::makeCorrectionVectors(vectorField& k) const { const vectorField unitDelta(delta()/mag(delta())); @@ -514,9 +532,10 @@ void Foam::faPatch::makeCorrectionVectors(vectorField& k) const } -const Foam::scalarField& Foam::faPatch::deltaCoeffs() const +const Foam::vectorField& Foam::faPatch::skewCorrectionVectors() const { - return boundaryMesh().mesh().deltaCoeffs().boundaryField()[index()]; + return + boundaryMesh().mesh().skewCorrectionVectors().boundaryField()[index()]; } diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H index 448b82d1e0..6bb87df10f 100644 --- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H +++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H @@ -413,14 +413,24 @@ public: //- Return patch weighting factors const scalarField& weights() const; + //- Make patch geodesic distance between P and N + virtual void makeLPN(scalarField&) const; + + //- Return patch geodesic distance between P and N + const scalarField& lPN() const; + //- Make patch edge - neighbour face distances virtual void makeDeltaCoeffs(scalarField&) const; - void makeCorrectionVectors(vectorField&) const; - //- Return patch edge - neighbour face distances const scalarField& deltaCoeffs() const; + //- Make non-orthogonality correction vectors + virtual void makeCorrectionVectors(vectorField&) const; + + //- Return patch skew-correction vectors + const vectorField& skewCorrectionVectors() const; + // Topological changes diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C index 5da45321ce..625a8a1471 100644 --- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C +++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C @@ -98,9 +98,17 @@ const Foam::edgeVectorField& Foam::edgeInterpolation::correctionVectors() const { if (orthogonal()) { - FatalErrorInFunction - << "cannot return correctionVectors; mesh is orthogonal" - << abort(FatalError); + return tmp::New + ( + IOobject + ( + "correctionVectors", + mesh().pointsInstance(), + mesh().thisDb() + ), + mesh(), + dimensionedVector(dimless, Zero) + ); } return (*correctionVectorsPtr_); @@ -123,9 +131,17 @@ Foam::edgeInterpolation::skewCorrectionVectors() const { if (!skew()) { - FatalErrorInFunction - << "cannot return skewCorrectionVectors; mesh is now skewed" - << abort(FatalError); + return tmp::New + ( + IOobject + ( + "skewCorrectionVectors", + mesh().pointsInstance(), + mesh().thisDb() + ), + mesh(), + dimensionedVector(dimless, Zero) + ); } return (*skewCorrectionVectorsPtr_); @@ -233,12 +249,10 @@ void Foam::edgeInterpolation::makeLPN() const forAll(lPN.boundaryField(), patchI) { - mesh().boundary()[patchI].makeDeltaCoeffs + mesh().boundary()[patchI].makeLPN ( lPN.boundaryFieldRef()[patchI] ); - - lPN.boundaryFieldRef()[patchI] = 1.0/lPN.boundaryField()[patchI]; } @@ -305,6 +319,7 @@ void Foam::edgeInterpolation::makeWeights() const // weight = (0,1] const scalar lPN = lPE + lEN; + if (mag(lPN) > SMALL) { weightingFactorsIn[edgeI] = lEN/lPN; @@ -498,31 +513,6 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const mesh().boundary()[patchI].makeCorrectionVectors(CorrVecsbf[patchI]); } - scalar NonOrthogCoeff = 0.0; - - if (owner.size() > 0) - { - scalarField sinAlpha(deltaCoeffs*mag(CorrVecs.internalField())); - sinAlpha.clamp_range(-1, 1); - - NonOrthogCoeff = max(Foam::asin(sinAlpha)*180.0/M_PI); - } - - reduce(NonOrthogCoeff, maxOp()); - - DebugInFunction - << "non-orthogonality coefficient = " << NonOrthogCoeff << " deg." - << endl; - - if (NonOrthogCoeff < 0.1) - { - orthogonal_ = true; - correctionVectorsPtr_.reset(nullptr); - } - else - { - orthogonal_ = false; - } DebugInFunction << "Finished constructing non-orthogonal correction vectors" diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoomWithThermalShell/system/finite-area/faSchemes b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoomWithThermalShell/system/finite-area/faSchemes index 8903be44f0..1810f4b3ce 100644 --- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoomWithThermalShell/system/finite-area/faSchemes +++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoomWithThermalShell/system/finite-area/faSchemes @@ -21,7 +21,7 @@ ddtSchemes gradSchemes { - default none; + default Gauss linear; grad(jouleHeatingSource:V_ceilingShell) Gauss linear; }