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.

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<symmTensor>& tf, const UList<symmTensor>& tf1)
+void inv(Field<symmTensor>& result, const UList<symmTensor>& tf1)
 {
-    if (tf.empty())
+    if (result.empty() || tf1.empty())
     {
         return;
     }
 
     // Attempting to identify 2-D cases
-    const scalar minThreshold = SMALL*magSqr(tf1[0]);
+    const scalar minThreshold = SMALL * magSqr(tf1[0]);
-    const Vector<bool> removeCmpts
-    (
-        magSqr(tf1[0].xx()) < minThreshold,
-        magSqr(tf1[0].yy()) < minThreshold,
-        magSqr(tf1[0].zz()) < minThreshold
-    );
 
-    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;
+        // Cannot use TFOR_ALL_F_OP_FUNC_F (additional operations)
-        if (removeCmpts.y()) adjust.yy() = 1;
-        if (removeCmpts.z()) adjust.zz() = 1;
 
-        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)
+            result[i] = Foam::inv(work);
-
+            result[i].subtractDiag(adjust);
-        tf -= 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<symmTensorField> inv(const tmp<symmTensorField>& tf)
     return tresult;
 }
 
+UNARY_FUNCTION(symmTensor, symmTensor, pinv)
+
 
 template<>
 tmp<Field<symmTensor>> transformFieldMask<symmTensor>
@@ -144,8 +159,6 @@ tmp<Field<symmTensor>> transformFieldMask<symmTensor>
     return tstf;
 }
 
-UNARY_FUNCTION(symmTensor, symmTensor, pinv)
-
 
 // * * * * * * * * * * * * * * * global operators  * * * * * * * * * * * * * //
 

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<tensor>& tf, const UList<tensor>& tf1)
+void inv(Field<tensor>& result, const UList<tensor>& tf1)
 {
-    if (tf.empty())
+    if (result.empty() || tf1.empty())
     {
         return;
     }
 
     // Attempting to identify 2-D cases
-    const scalar minThreshold = SMALL*magSqr(tf1[0]);
+    const scalar minThreshold = SMALL * magSqr(tf1[0]);
-    const Vector<bool> removeCmpts
-    (
-        magSqr(tf1[0].xx()) < minThreshold,
-        magSqr(tf1[0].yy()) < minThreshold,
-        magSqr(tf1[0].zz()) < minThreshold
-    );
 
-    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;
+        // Cannot use TFOR_ALL_F_OP_FUNC_F (additional operations)
-        if (removeCmpts.y()) adjust.yy() = 1;
-        if (removeCmpts.z()) adjust.zz() = 1;
 
-        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)
+            result[i] = Foam::inv(work);
-
+            result[i].subtractDiag(adjust);
-        tf -= 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]);
+        }
     }
 }
 

src/OpenFOAM/primitives/SymmTensor/SymmTensor.H

@@ -251,8 +251,9 @@ public:
         //- Return 2D cofactor matrix by excluding given direction
         inline SymmTensor<Cmpt> cof2D(const direction excludeCmpt) const;
 
-        //- Return inverse
+        //- Return inverse, with optional (ad hoc) failsafe handling
-        inline SymmTensor<Cmpt> inv() const;
+        //- of 2D tensors
+        inline SymmTensor<Cmpt> inv(const bool failsafe=false) const;
 
         //- Return inverse of 2D tensor (by excluding given direction)
         inline SymmTensor<Cmpt> inv2D(const direction excludeCmpt) const;

src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H

@@ -451,9 +451,57 @@ inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st)
 
 
 template<class Cmpt>
-inline SymmTensor<Cmpt> SymmTensor<Cmpt>::inv() const
+inline SymmTensor<Cmpt> SymmTensor<Cmpt>::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<Cmpt> work(*this);
+
+            if (small_xx) { work.xx() += pTraits<Cmpt>::one; }
+            if (small_yy) { work.yy() += pTraits<Cmpt>::one; }
+            if (small_zz) { work.zz() += pTraits<Cmpt>::one; }
+
+            work = Foam::inv(work, work.det());
+
+            if (small_xx) { work.xx() -= pTraits<Cmpt>::one; }
+            if (small_yy) { work.yy() -= pTraits<Cmpt>::one; }
+            if (small_zz) { work.zz() -= pTraits<Cmpt>::one; }
+
+            return work;
+        }
+        else if (mag(detval) < VSMALL)
+        {
+            // Really appears to be zero - leave untouched?
+            return SymmTensor<Cmpt>(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;
 }
 
 

src/OpenFOAM/primitives/Tensor/Tensor.H

@@ -293,8 +293,9 @@ public:
         //- Return 2D cofactor matrix by excluding given direction
         inline Tensor<Cmpt> cof2D(const direction excludeCmpt) const;
 
-        //- Return inverse
+        //- Return inverse, with optional (ad hoc) failsafe handling
-        inline Tensor<Cmpt> inv() const;
+        //- of 2D tensors
+        inline Tensor<Cmpt> inv(const bool failsafe=false) const;
 
         //- Return inverse of 2D tensor (by excluding given direction)
         inline Tensor<Cmpt> inv2D(const direction excludeCmpt) const;

src/OpenFOAM/primitives/Tensor/TensorI.H

@@ -762,9 +762,58 @@ inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t)
 
 
 template<class Cmpt>
-inline Tensor<Cmpt> Tensor<Cmpt>::inv() const
+inline Tensor<Cmpt> Tensor<Cmpt>::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<Cmpt> work(*this);
+
+            if (small_xx) { work.xx() += pTraits<Cmpt>::one; }
+            if (small_yy) { work.yy() += pTraits<Cmpt>::one; }
+            if (small_zz) { work.zz() += pTraits<Cmpt>::one; }
+
+            work = Foam::inv(work, work.det());
+
+            if (small_xx) { work.xx() -= pTraits<Cmpt>::one; }
+            if (small_yy) { work.yy() -= pTraits<Cmpt>::one; }
+            if (small_zz) { work.zz() -= pTraits<Cmpt>::one; }
+
+            return work;
+        }
+        else if (mag(detval) < VSMALL)
+        {
+            // Really appears to be zero - leave untouched?
+            return Tensor<Cmpt>(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;
 }
 
 

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
+    // Invert the dd tensor.
-    const symmTensorField invDd(inv(dd));
+
+    // 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