From 66cae2b9d124caad78ea6d2321d8e85343aec67b Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Wed, 22 Mar 2023 22:58:31 +0100 Subject: [PATCH] ENH: robuster handling inv() of singular tensor for finite-area LSQ (#2724) - with the current handling of small edges (finite-area), the LSQ vectors can result in singular/2D tensors. However, the regular 2D handling in field inv() only detects based on the first element. Provide a 'failsafe' inv() method for symmTensor and tensor that follows a similar logic for avoiding zero determinates, but it is applied on a per element basis, instead of deciding based on the first field element. The symmTensor::inv(bool) and tensor::inv(bool) methods have a fairly modest additional overhead. - unroll the field inv() function to avoid creating an intermediate field. Reduce the number of operations when adjusting/re-adjusting the diagonal. --- .../Fields/symmTensorField/symmTensorField.C | 59 +++++++++++-------- .../fields/Fields/tensorField/tensorField.C | 53 ++++++++++------- .../primitives/SymmTensor/SymmTensor.H | 5 +- .../primitives/SymmTensor/SymmTensorI.H | 52 +++++++++++++++- src/OpenFOAM/primitives/Tensor/Tensor.H | 5 +- src/OpenFOAM/primitives/Tensor/TensorI.H | 53 ++++++++++++++++- .../leastSquaresFaVectors.C | 28 +++++++-- 7 files changed, 199 insertions(+), 56 deletions(-) diff --git a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C index d71bafa5e5..1be5183d4c 100644 --- a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C +++ b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019-2022 OpenCFD Ltd. + Copyright (C) 2019-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -51,41 +51,54 @@ UNARY_FUNCTION(symmTensor, symmTensor, dev2) UNARY_FUNCTION(scalar, symmTensor, det) UNARY_FUNCTION(symmTensor, symmTensor, cof) -void inv(Field& tf, const UList& tf1) +void inv(Field& result, const UList& tf1) { - if (tf.empty()) + if (result.empty() || tf1.empty()) { return; } // Attempting to identify 2-D cases - const scalar minThreshold = SMALL*magSqr(tf1[0]); - const Vector removeCmpts - ( - magSqr(tf1[0].xx()) < minThreshold, - magSqr(tf1[0].yy()) < minThreshold, - magSqr(tf1[0].zz()) < minThreshold - ); + const scalar minThreshold = SMALL * magSqr(tf1[0]); - if (removeCmpts.x() || removeCmpts.y() || removeCmpts.z()) + const bool small_xx = (magSqr(tf1[0].xx()) < minThreshold); + const bool small_yy = (magSqr(tf1[0].yy()) < minThreshold); + const bool small_zz = (magSqr(tf1[0].zz()) < minThreshold); + + if (small_xx || small_yy || small_zz) { - symmTensor adjust(Zero); + const vector adjust + ( + (small_xx ? 1 : 0), + (small_yy ? 1 : 0), + (small_zz ? 1 : 0) + ); - if (removeCmpts.x()) adjust.xx() = 1; - if (removeCmpts.y()) adjust.yy() = 1; - if (removeCmpts.z()) adjust.zz() = 1; + // Cannot use TFOR_ALL_F_OP_FUNC_F (additional operations) - symmTensorField tf1Plus(tf1); + const label loopLen = (result).size(); - tf1Plus += adjust; + /* pragmas... */ + for (label i = 0; i < loopLen; ++i) + { + symmTensor work(tf1[i]); + work.addDiag(adjust); - TFOR_ALL_F_OP_FUNC_F(symmTensor, tf, =, inv, symmTensor, tf1Plus) - - tf -= adjust; + result[i] = Foam::inv(work); + result[i].subtractDiag(adjust); + } } else { - TFOR_ALL_F_OP_FUNC_F(symmTensor, tf, =, inv, symmTensor, tf1) + // Same as TFOR_ALL_F_OP_FUNC_F + + const label loopLen = (result).size(); + + /* pragmas... */ + for (label i = 0; i < loopLen; ++i) + { + result[i] = Foam::inv(tf1[i]); + } } } @@ -104,6 +117,8 @@ tmp inv(const tmp& tf) return tresult; } +UNARY_FUNCTION(symmTensor, symmTensor, pinv) + template<> tmp> transformFieldMask @@ -144,8 +159,6 @@ tmp> transformFieldMask return tstf; } -UNARY_FUNCTION(symmTensor, symmTensor, pinv) - // * * * * * * * * * * * * * * * global operators * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C index 43eb0b8493..c4f76fc978 100644 --- a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C +++ b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C @@ -50,41 +50,54 @@ UNARY_FUNCTION(tensor, tensor, dev2) UNARY_FUNCTION(scalar, tensor, det) UNARY_FUNCTION(tensor, tensor, cof) -void inv(Field& tf, const UList& tf1) +void inv(Field& result, const UList& tf1) { - if (tf.empty()) + if (result.empty() || tf1.empty()) { return; } // Attempting to identify 2-D cases - const scalar minThreshold = SMALL*magSqr(tf1[0]); - const Vector removeCmpts - ( - magSqr(tf1[0].xx()) < minThreshold, - magSqr(tf1[0].yy()) < minThreshold, - magSqr(tf1[0].zz()) < minThreshold - ); + const scalar minThreshold = SMALL * magSqr(tf1[0]); - if (removeCmpts.x() || removeCmpts.y() || removeCmpts.z()) + const bool small_xx = (magSqr(tf1[0].xx()) < minThreshold); + const bool small_yy = (magSqr(tf1[0].yy()) < minThreshold); + const bool small_zz = (magSqr(tf1[0].zz()) < minThreshold); + + if (small_xx || small_yy || small_zz) { - tensor adjust(Zero); + const vector adjust + ( + (small_xx ? 1 : 0), + (small_yy ? 1 : 0), + (small_zz ? 1 : 0) + ); - if (removeCmpts.x()) adjust.xx() = 1; - if (removeCmpts.y()) adjust.yy() = 1; - if (removeCmpts.z()) adjust.zz() = 1; + // Cannot use TFOR_ALL_F_OP_FUNC_F (additional operations) - tensorField tf1Plus(tf1); + const label loopLen = (result).size(); - tf1Plus += adjust; + /* pragmas... */ + for (label i = 0; i < loopLen; ++i) + { + tensor work(tf1[i]); + work.addDiag(adjust); - TFOR_ALL_F_OP_FUNC_F(tensor, tf, =, inv, tensor, tf1Plus) - - tf -= adjust; + result[i] = Foam::inv(work); + result[i].subtractDiag(adjust); + } } else { - TFOR_ALL_F_OP_FUNC_F(tensor, tf, =, inv, tensor, tf1) + // Same as TFOR_ALL_F_OP_FUNC_F + + const label loopLen = (result).size(); + + /* pragmas... */ + for (label i = 0; i < loopLen; ++i) + { + result[i] = Foam::inv(tf1[i]); + } } } diff --git a/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H b/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H index d3f3f5f88f..e365081023 100644 --- a/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H +++ b/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H @@ -251,8 +251,9 @@ public: //- Return 2D cofactor matrix by excluding given direction inline SymmTensor cof2D(const direction excludeCmpt) const; - //- Return inverse - inline SymmTensor inv() const; + //- Return inverse, with optional (ad hoc) failsafe handling + //- of 2D tensors + inline SymmTensor inv(const bool failsafe=false) const; //- Return inverse of 2D tensor (by excluding given direction) inline SymmTensor inv2D(const direction excludeCmpt) const; diff --git a/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H b/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H index 56ad9eb4c6..47c10fa13b 100644 --- a/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H +++ b/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H @@ -451,9 +451,57 @@ inline SymmTensor inv(const SymmTensor& st) template -inline SymmTensor SymmTensor::inv() const +inline SymmTensor SymmTensor::inv(const bool failsafe) const { - return Foam::inv(*this, this->det()); + const Cmpt detval = this->det(); + + if (failsafe) + { + // Attempt to identify and handle 2-D cases. + // - use diagSqr instead of magSqr for fewer operations + const scalar magSqr_xx = Foam::magSqr(xx()); + const scalar magSqr_yy = Foam::magSqr(yy()); + const scalar magSqr_zz = Foam::magSqr(zz()); + + const scalar threshold = SMALL * (magSqr_xx + magSqr_yy + magSqr_zz); + + const bool small_xx = (magSqr_xx < threshold); + const bool small_yy = (magSqr_yy < threshold); + const bool small_zz = (magSqr_zz < threshold); + + if (small_xx || small_yy || small_zz) + { + SymmTensor work(*this); + + if (small_xx) { work.xx() += pTraits::one; } + if (small_yy) { work.yy() += pTraits::one; } + if (small_zz) { work.zz() += pTraits::one; } + + work = Foam::inv(work, work.det()); + + if (small_xx) { work.xx() -= pTraits::one; } + if (small_yy) { work.yy() -= pTraits::one; } + if (small_zz) { work.zz() -= pTraits::one; } + + return work; + } + else if (mag(detval) < VSMALL) + { + // Really appears to be zero - leave untouched? + return SymmTensor(Zero); + } + } + #ifdef FULLDEBUG + else if (mag(detval) < VSMALL) + { + FatalErrorInFunction + << "Tensor not properly invertible, determinant:" + << detval << " tensor:" << *this << nl + << abort(FatalError); + } + #endif + + return this->cof().T()/detval; } diff --git a/src/OpenFOAM/primitives/Tensor/Tensor.H b/src/OpenFOAM/primitives/Tensor/Tensor.H index eac4435f94..7391d06865 100644 --- a/src/OpenFOAM/primitives/Tensor/Tensor.H +++ b/src/OpenFOAM/primitives/Tensor/Tensor.H @@ -293,8 +293,9 @@ public: //- Return 2D cofactor matrix by excluding given direction inline Tensor cof2D(const direction excludeCmpt) const; - //- Return inverse - inline Tensor inv() const; + //- Return inverse, with optional (ad hoc) failsafe handling + //- of 2D tensors + inline Tensor inv(const bool failsafe=false) const; //- Return inverse of 2D tensor (by excluding given direction) inline Tensor inv2D(const direction excludeCmpt) const; diff --git a/src/OpenFOAM/primitives/Tensor/TensorI.H b/src/OpenFOAM/primitives/Tensor/TensorI.H index a0254112bb..ec6997262a 100644 --- a/src/OpenFOAM/primitives/Tensor/TensorI.H +++ b/src/OpenFOAM/primitives/Tensor/TensorI.H @@ -762,9 +762,58 @@ inline Tensor inv(const Tensor& t) template -inline Tensor Tensor::inv() const +inline Tensor Tensor::inv(const bool failsafe) const { - return Foam::inv(*this, this->det()); + const Cmpt detval = this->det(); + + if (failsafe) + { + // Attempt to identify and handle 2-D cases + // - use diagSqr instead of magSqr for fewer operations + + const scalar magSqr_xx = Foam::magSqr(xx()); + const scalar magSqr_yy = Foam::magSqr(yy()); + const scalar magSqr_zz = Foam::magSqr(zz()); + + const scalar threshold = SMALL * (magSqr_xx + magSqr_yy + magSqr_zz); + + const bool small_xx = (magSqr_xx < threshold); + const bool small_yy = (magSqr_yy < threshold); + const bool small_zz = (magSqr_zz < threshold); + + if (small_xx || small_yy || small_zz) + { + Tensor work(*this); + + if (small_xx) { work.xx() += pTraits::one; } + if (small_yy) { work.yy() += pTraits::one; } + if (small_zz) { work.zz() += pTraits::one; } + + work = Foam::inv(work, work.det()); + + if (small_xx) { work.xx() -= pTraits::one; } + if (small_yy) { work.yy() -= pTraits::one; } + if (small_zz) { work.zz() -= pTraits::one; } + + return work; + } + else if (mag(detval) < VSMALL) + { + // Really appears to be zero - leave untouched? + return Tensor(Zero); + } + } + #ifdef FULLDEBUG + else if (mag(detval) < VSMALL) + { + FatalErrorInFunction + << "Tensor not properly invertible, determinant:" + << detval << " tensor:" << *this << nl + << abort(FatalError); + } + #endif + + return this->cof().T()/detval; } diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C index ebcd3d6d64..f2bd352bbd 100644 --- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C +++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2016-2017 Wikki Ltd - Copyright (C) 2020-2022 OpenCFD Ltd. + Copyright (C) 2020-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -81,7 +81,7 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const mesh(), dimensionedVector(dimless/dimLength, Zero) ); - edgeVectorField& lsP = *pVectorsPtr_; + auto& lsP = *pVectorsPtr_; nVectorsPtr_ = new edgeVectorField ( @@ -97,7 +97,7 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const mesh(), dimensionedVector(dimless/dimLength, Zero) ); - edgeVectorField& lsN = *nVectorsPtr_; + auto& lsN = *nVectorsPtr_; // Set local references to mesh data const labelUList& owner = mesh().owner(); @@ -148,8 +148,26 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const } - // Invert the dd tensor - const symmTensorField invDd(inv(dd)); + // Invert the dd tensor. + + // Cannot rely on the usual field inv() since that only uses the + // first element to guess if the remaining tensors are 2D (or + // singular). We, however, can have a mixture (eg, skipping over + // zero-length edges can yield a zero). Instead use the + // 'failsafe' inv() that checks each tensor (#2724) + + // Fragile: const symmTensorField invDd(inv(dd)); + + symmTensorField invDd(dd.size()); + { + const label loopLen = dd.size(); + + /* pragmas... */ + for (label i = 0; i < loopLen; ++i) + { + invDd[i] = dd[i].inv(true); // With 'failsafe' handling + } + } // Revisit all faces and calculate the lsP and lsN vectors