diff --git a/applications/test/SymmTensor/Test-SymmTensor.C b/applications/test/SymmTensor/Test-SymmTensor.C
index 2d1aaa970a..35da47cd8a 100644
--- a/applications/test/SymmTensor/Test-SymmTensor.C
+++ b/applications/test/SymmTensor/Test-SymmTensor.C
@@ -62,9 +62,24 @@ unsigned nFail_ = 0;
 // Create a random symmTensor
 symmTensor makeRandomContainer(Random& rnd)
 {
-    symmTensor A(Zero);
-    std::generate(A.begin(), A.end(), [&]{ return rnd.GaussNormal<scalar>(); });
-    return A;
+    symmTensor T(Zero);
+    std::generate(T.begin(), T.end(), [&]{ return rnd.GaussNormal<scalar>(); });
+    return T;
+}
+
+
+// Create a symmTensor based on a given value
+template<class Type>
+typename std::enable_if
+<
+    std::is_same<floatScalar, Type>::value ||
+    std::is_same<doubleScalar, Type>::value,
+    symmTensor
+>::type makeContainer(const Type val)
+{
+    symmTensor T(Zero);
+    std::fill(T.begin(), T.end(), val);
+    return T;
 }
 
 
@@ -83,11 +98,11 @@ typename std::enable_if
     const word& msg,
     const Type& x,
     const Type& y,
-    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
-    const scalar absTol = 0      //<! useful for cmps near zero
+    const scalar absTol = 0,    //<! useful for cmps near zero
+    const scalar relTol = 1e-8  //<! are values the same within 8 decimals
 )
 {
-    Info<< msg << x << endl;
+    Info<< msg << x << "?=" << y << endl;
 
     unsigned nFail = 0;
 
@@ -121,11 +136,11 @@ typename std::enable_if
     const word& msg,
     const Type& x,
     const Type& y,
-    const scalar relTol = 1e-8,
-    const scalar absTol = 0
+    const scalar absTol = 0,
+    const scalar relTol = 1e-8
 )
 {
-    Info<< msg << x << endl;
+    Info<< msg << x << "?=" << y << endl;
 
     unsigned nFail = 0;
 
@@ -338,7 +353,6 @@ void test_global_funcs(Type)
                        Type(1),  Type(0),
                                  Type(-0.05555556)
         ),
-        1e-8,
         1e-8
     );
     cmp
@@ -351,7 +365,6 @@ void test_global_funcs(Type)
                        Type(1),  Type(0),
                                  Type(-0.05555556)
         ),
-        1e-8,
         1e-8
     );
     cmp("  First invariant = ", invariantI(sT), Type(-3));
@@ -514,21 +527,21 @@ void test_global_opers(Type)
 
 // Return false if given eigenvalues fail to satisy eigenvalue relations
 // Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307)
-void test_eigenvalues(const symmTensor& sT, const vector& EVals)
+void test_eigenvalues(const symmTensor& T, const vector& EVals)
 {
     {
-        const scalar determinant = det(sT);
+        const scalar determinant = det(T);
         const scalar EValsProd = EVals.x()*EVals.y()*EVals.z();
-        cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant);
+        cmp("# Product of eigenvalues = det(T):", EValsProd, determinant, 1e-6);
     }
     {
-        const scalar trace = tr(sT);
+        const scalar trace = tr(T);
         scalar EValsSum = 0.0;
         for (const auto& val : EVals)
         {
             EValsSum += val;
         }
-        cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace);
+        cmp("# Sum of eigenvalues = trace(T):", EValsSum, trace);
     }
 }
 
@@ -537,7 +550,7 @@ void test_eigenvalues(const symmTensor& sT, const vector& EVals)
 // fails to satisfy the characteristic equation
 void test_characteristic_equation
 (
-    const symmTensor& sT,
+    const symmTensor& T,
     const vector& EVals,
     const tensor& EVecs
 )
@@ -547,39 +560,39 @@ void test_characteristic_equation
     {
         Info<< "EVal = " << EVals[dir] << nl
             << "EVec = " << EVecs.row(dir) << endl;
-        const vector leftSide(sT & EVecs.row(dir));
+        const vector leftSide(T & EVecs.row(dir));
         const vector rightSide(EVals[dir]*EVecs.row(dir));
         const vector X(leftSide - rightSide);
 
         for (const auto x : X)
         {
-            cmp("  (sT & EVec - EVal*EVec) = 0:", x, 0.0, 1e-8, 1e-6);
+            cmp("  (T & EVec - EVal*EVec) = 0:", x, 0.0, 1e-5);
         }
     }
 }
 
 
 // Return false if the eigen functions fail to satisfy relations
-void test_eigen_funcs(const symmTensor& sT)
+void test_eigen_funcs(const symmTensor& T)
 {
     Info<< "# Operand:" << nl
-        << "  symmTensor = " << sT << nl;
+        << "  symmTensor = " << T << nl;
 
 
     Info<< "# Return eigenvalues of a given symmTensor:" << nl;
-    const vector EVals(eigenValues(sT));
+    const vector EVals(eigenValues(T));
     Info<< EVals << endl;
-    test_eigenvalues(sT, EVals);
+    test_eigenvalues(T, EVals);
 
     Info<< "# Return eigenvectors of a given symmTensor corresponding to"
         << " given eigenvalues:" << nl;
-    const tensor EVecs0(eigenVectors(sT, EVals));
+    const tensor EVecs0(eigenVectors(T, EVals));
     Info<< EVecs0 << endl;
-    test_characteristic_equation(sT, EVals, EVecs0);
+    test_characteristic_equation(T, EVals, EVecs0);
 
     Info<< "# Return eigenvectors of a given symmTensor by computing"
         << " the eigenvalues of the symmTensor in the background:" << nl;
-    const tensor EVecs1(eigenVectors(sT));
+    const tensor EVecs1(eigenVectors(T));
     Info<< EVecs1 << endl;
 }
 
@@ -629,7 +642,7 @@ int main()
     run_tests(types, typeID);
 
 
-    Info<< nl << "    ## Test SymmTensor<scalar> eigen functions: ##" << nl;
+    Info<< nl << "    ## Test symmTensor eigen functions: ##" << nl;
     const label numberOfTests = 10000;
     Random rndGen(1234);
 
@@ -640,14 +653,13 @@ int main()
     }
 
     {
-        Info<< nl << "    ## Test eigen functions by a zero tensor: ##"<< nl;
+        Info<< nl << "    ## A zero symmTensor: ##"<< nl;
         const symmTensor zeroT(Zero);
         test_eigen_funcs(zeroT);
     }
     {
         Info<< nl
-            << "    ## Test eigen functions by a tensor consisting of"
-            << " 2 repeated eigenvalues: ##"
+            << "    ## A symmTensor with 2 repeated eigenvalues: ##"
             << nl;
         const symmTensor T
         (
@@ -659,8 +671,7 @@ int main()
     }
     {
         Info<< nl
-            << "    ## Test eigen functions by a tensor consisting of"
-            << " 3 repeated eigenvalues: ##"
+            << "    ## A symmTensor with 3 repeated eigenvalues: ##"
             << nl;
         const symmTensor T
         (
@@ -671,22 +682,7 @@ int main()
         test_eigen_funcs(T);
     }
     {
-        Info<< nl
-            << "    ## Test eigen functions by a tensor consisting of"
-            << " SMALL off-diagonal elements: ##"
-            << nl;
-
-        for (label i = 0; i < numberOfTests; ++i)
-        {
-            symmTensor T(makeRandomContainer(rndGen));
-            T.xy() = SMALL*rndGen.GaussNormal<scalar>();
-            T.xz() = SMALL*rndGen.GaussNormal<scalar>();
-            T.yz() = SMALL*rndGen.GaussNormal<scalar>();
-            test_eigen_funcs(T);
-        }
-    }
-    {
-        Info<< nl << "    ## Test eigen functions by a stiff tensor: ##" << nl;
+        Info<< nl << "    ## A stiff symmTensor: ##" << nl;
         const symmTensor stiff
         (
             pow(10.0, 10), pow(10.0, 8),  pow(10.0, -8),
@@ -695,6 +691,165 @@ int main()
         );
         test_eigen_funcs(stiff);
     }
+    {
+        Info<< nl
+            << "    ## Random symmTensors with tiny off-diag elements: ##"
+            << nl;
+
+        const List<scalar> epsilons
+        ({
+            0, SMALL, Foam::sqrt(SMALL), sqr(SMALL), Foam::cbrt(SMALL),
+            -SMALL, -Foam::sqrt(SMALL), -sqr(SMALL), -Foam::cbrt(SMALL)
+        });
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            for (const auto& eps : epsilons)
+            {
+                {
+                    symmTensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    symmTensor T(makeRandomContainer(rndGen));
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    symmTensor T(makeRandomContainer(rndGen));
+                    T.yz() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    symmTensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    symmTensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    T.yz() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    symmTensor T(makeRandomContainer(rndGen));
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yz() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    symmTensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yz() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    symmTensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps;
+                    T.xz() = eps;
+                    T.yz() = eps;
+                    test_eigen_funcs(T);
+                }
+                {
+                    symmTensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps;
+                    T.xz() = eps;
+                    T.yz() = eps;
+                    T.zz() = eps;
+                    test_eigen_funcs(T);
+                }
+                {
+                    symmTensor T(makeRandomContainer(rndGen));
+                    T.xy() = 0;
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yz() = 0;
+                    test_eigen_funcs(T);
+                }
+            }
+        }
+    }
+    #if 0
+    // Numerical diagonalisation of 2x2 or 3x3 matrices with analytic methods
+    // are, like the methods currently being used in OpenFOAM, inherently error
+    // prone. Despite its speed, the analytic methods may becomes inaccurate or
+    // may even fail completely if the matrix entries differ greatly in
+    // magnitude, particularly with large off-diagonal elements.
+    // The remedy is to use iterative or hybrid analytic/iterative methods
+    // such as published here (for 3x3/2x2 matrices):
+    // (Kopp, 2008) arXiv.org: physics/0610206
+    // mpi-hd.mpg.de/personalhomes/globes/3x3/index.html
+    {
+        Info<< nl << "    ## symmTensors consisting machine epsilons: ##" << nl;
+        Info<< "    # floatScalar" << nl;
+        const List<floatScalar> floatEpsilons
+        ({
+            floatScalarGREAT, floatScalarVGREAT, floatScalarROOTVGREAT,
+            floatScalarSMALL, floatScalarVSMALL, floatScalarROOTVSMALL,
+            Foam::sqrt(floatScalarSMALL), 0
+        });
+
+        for (const auto& eps : floatEpsilons)
+        {
+            const symmTensor T(makeContainer(eps));
+            test_eigen_funcs(T);
+        }
+
+        Info<< "    # doubleScalar" << nl;
+        const List<doubleScalar> doubleEpsilons
+        ({
+            doubleScalarGREAT, doubleScalarROOTVGREAT, // doubleVGREAT fails
+            doubleScalarSMALL, doubleScalarVSMALL, doubleScalarROOTVSMALL,
+            Foam::sqrt(doubleScalarSMALL), 0
+        });
+
+        for (const auto& eps : doubleEpsilons)
+        {
+            const symmTensor T(makeContainer(eps));
+            test_eigen_funcs(T);
+        }
+    }
+    {
+        Info<< nl
+            << "    ## Random symmTensors with machine eps off-diag elmes: ##"
+            << nl;
+
+        const List<floatScalar> floatEpsilons
+        ({
+            floatScalarGREAT, floatScalarVGREAT, floatScalarROOTVGREAT,
+            floatScalarSMALL, floatScalarVSMALL, floatScalarROOTVSMALL
+        });
+
+        const List<doubleScalar> doubleEpsilons
+        ({
+            doubleScalarGREAT, doubleScalarVGREAT, doubleScalarROOTVGREAT,
+            doubleScalarSMALL, doubleScalarVSMALL, doubleScalarROOTVSMALL
+        });
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            symmTensor T(makeRandomContainer(rndGen));
+
+            for (const auto& eps : floatEpsilons)
+            {
+                T.xy() = eps;
+                T.xz() = eps;
+                T.yz() = eps;
+                test_eigen_funcs(T);
+            }
+
+            for (const auto& eps : doubleEpsilons)
+            {
+                T.xy() = eps;
+                T.xz() = eps;
+                T.yz() = eps;
+                test_eigen_funcs(T);
+            }
+        }
+    }
+    #endif
 
 
     if (nFail_)
diff --git a/applications/test/SymmTensor2D/Test-SymmTensor2D.C b/applications/test/SymmTensor2D/Test-SymmTensor2D.C
index 82bbc77745..8b0485132c 100644
--- a/applications/test/SymmTensor2D/Test-SymmTensor2D.C
+++ b/applications/test/SymmTensor2D/Test-SymmTensor2D.C
@@ -83,11 +83,11 @@ typename std::enable_if
     const word& msg,
     const Type& x,
     const Type& y,
-    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
-    const scalar absTol = 0      //<! useful for cmps near zero
+    const scalar absTol = 0,     //<! useful for cmps near zero
+    const scalar relTol = 1e-8   //<! are values the same within 8 decimals
 )
 {
-    Info<< msg << x << endl;
+    Info<< msg << x << "?=" << y << endl;
 
     unsigned nFail = 0;
 
@@ -121,11 +121,11 @@ typename std::enable_if
     const word& msg,
     const Type& x,
     const Type& y,
-    const scalar relTol = 1e-8,
-    const scalar absTol = 0
+    const scalar absTol = 0,
+    const scalar relTol = 1e-8
 )
 {
-    Info<< msg << x << endl;
+    Info<< msg << x << "?=" << y << endl;
 
     unsigned nFail = 0;
 
@@ -311,8 +311,7 @@ void test_global_funcs(Type)
             Type(0.42857143), Type(0.28571429),
                               Type(-0.14285714)
         ),
-        1e-8,
-        1e-6
+        1e-7
     );
     cmp
     (
@@ -323,8 +322,7 @@ void test_global_funcs(Type)
             Type(0.42857143), Type(0.28571429),
                               Type(-0.14285714)
         ),
-        1e-8,
-        1e-6
+        1e-7
     );
     cmp("  First invariant = ", invariantI(sT), Type(-2));
     cmp("  Second invariant = ", invariantII(sT), Type(-7));
@@ -481,10 +479,9 @@ void test_global_opers(Type)
 void test_eigenvalues(const symmTensor2D& T, const vector2D& EVals)
 {
     {
-        Info<< "# Product of eigenvalues = det(T):" << nl;
         const scalar determinant = det(T);
         const scalar EValsProd = EVals.x()*EVals.y();
-        cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant);
+        cmp("# Product of eigenvalues = det(T):", EValsProd, determinant, 1e-8);
     }
 
     {
@@ -494,7 +491,7 @@ void test_eigenvalues(const symmTensor2D& T, const vector2D& EVals)
         {
             EValsSum += val;
         }
-        cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace);
+        cmp("# Sum of eigenvalues = trace(T):", EValsSum, trace);
     }
 }
 
@@ -519,7 +516,7 @@ void test_characteristic_equation
 
         for (const auto x : X)
         {
-            cmp("  (sT & EVec - EVal*EVec) = 0:", x, 0.0, 1e-8, 1e-6);
+            cmp("  (sT & EVec - EVal*EVec) = 0:", x, 0.0, 1e-5);
         }
     }
 }
@@ -605,33 +602,31 @@ int main(int argc, char *argv[])
         test_eigen_funcs(T);
     }
 
-    Info<< nl << "    ## Test symmTensor2D eigen functions"
-        << " with T.xy = VSMALL: ##" << nl;
-    for (label i = 0; i < numberOfTests; ++i)
     {
-        symmTensor2D T(makeRandomContainer(rndGen));
-        T.xy() = VSMALL;
-        test_eigen_funcs(T);
+        Info<< nl << "    ## A symmTensor2D with T.xy = VSMALL" << nl;
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            symmTensor2D T(makeRandomContainer(rndGen));
+            T.xy() = VSMALL;
+            test_eigen_funcs(T);
+        }
     }
-
-    Info<< nl << "    ## Test symmTensor2D eigen functions"
-        << " with T.xx = T.yy: ##" << nl;
-    for (label i = 0; i < numberOfTests; ++i)
     {
-        symmTensor2D T(makeRandomContainer(rndGen));
-        T.xx() = T.yy();
-        test_eigen_funcs(T);
+        Info<< nl << "    ## A symmTensor2D with T.xx = T.yy: ##" << nl;
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            symmTensor2D T(makeRandomContainer(rndGen));
+            T.xx() = T.yy();
+            test_eigen_funcs(T);
+        }
     }
-
     {
-        Info<< nl
-            << "    ## Test eigen functions by a zero symmTensor2D: ##"<< nl;
+        Info<< nl << "    ## A zero symmTensor2D: ##"<< nl;
         const symmTensor2D zeroT(Zero);
         test_eigen_funcs(zeroT);
     }
-
     {
-        Info<< nl << "    ## Test eigen functions by a stiff tensor2D: ##"<< nl;
+        Info<< nl << "    ## A stiff symmTensor2D: ##"<< nl;
         const symmTensor2D stiff
         (
             pow(10.0, 10), pow(10.0, -8),
@@ -639,6 +634,34 @@ int main(int argc, char *argv[])
         );
         test_eigen_funcs(stiff);
     }
+    {
+        Info<< nl
+            << "    ## Random symmTensor2D with tiny off-diag elements: ##"
+            << nl;
+
+        const List<scalar> epsilons
+        ({
+            0, SMALL, Foam::sqrt(SMALL), sqr(SMALL), Foam::cbrt(SMALL),
+            -SMALL, -Foam::sqrt(SMALL), -sqr(SMALL), -Foam::cbrt(SMALL)
+        });
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            for (const auto& eps : epsilons)
+            {
+                {
+                    symmTensor2D T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    symmTensor2D T(makeRandomContainer(rndGen));
+                    T.xy() = eps;
+                    test_eigen_funcs(T);
+                }
+            }
+        }
+    }
 
 
     if (nFail_)
diff --git a/applications/test/Tensor/Test-Tensor.C b/applications/test/Tensor/Test-Tensor.C
index 18f8d2b51c..d9142c467d 100644
--- a/applications/test/Tensor/Test-Tensor.C
+++ b/applications/test/Tensor/Test-Tensor.C
@@ -84,11 +84,11 @@ typename std::enable_if
     const word& msg,
     const Type& x,
     const Type& y,
-    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
-    const scalar absTol = 0      //<! useful for cmps near zero
+    const scalar absTol = 0,     //<! useful for cmps near zero
+    const scalar relTol = 1e-8   //<! are values the same within 8 decimals
 )
 {
-    Info<< msg << x << endl;
+    Info<< msg << x << "?=" << y << endl;
 
     unsigned nFail = 0;
 
@@ -122,11 +122,11 @@ typename std::enable_if
     const word& msg,
     const Type& x,
     const Type& y,
-    const scalar relTol = 1e-8,
-    const scalar absTol = 0
+    const scalar absTol = 0,
+    const scalar relTol = 1e-8
 )
 {
-    Info<< msg << x << endl;
+    Info<< msg << x << "?=" << y << endl;
 
     unsigned nFail = 0;
 
@@ -601,8 +601,7 @@ void test_global_funcs(Type)
             Type(4),          Type(6.66666667), Type(-6),
             Type(7),          Type(8),          Type(-7.33333333)
         ),
-        1e-6,
-        1e-6
+        1e-7
     );
     cmp("  Two-third deviatoric part = ", dev2(T), T - 2*sph(T));
     cmp("  Determinant = ", det(T), Type(-6.000000000000005));
@@ -627,7 +626,6 @@ void test_global_funcs(Type)
             Type(1),     Type(-5),           Type(3),
             Type(0.5),   Type(-3.66666667),  Type(2.16666667)
         ),
-        1e-8,
         1e-8
     );
     cmp
@@ -640,7 +638,6 @@ void test_global_funcs(Type)
             Type(1),     Type(-5),           Type(3),
             Type(0.5),   Type(-3.66666667),  Type(2.16666667)
         ),
-        1e-8,
         1e-8
     );
     cmp
@@ -653,7 +650,6 @@ void test_global_funcs(Type)
             Type(1),     Type(-5),           Type(3),
             Type(0.5),   Type(-3.66666667),  Type(2.16666667)
         ),
-        1e-8,
         1e-8
     );
     cmp("  First invariant = ", invariantI(T), Type(-5));
@@ -917,7 +913,7 @@ void test_eigenvalues
         // In case of complex EVals, the production is effectively scalar
         // due to the (complex*complex conjugate) results in zero imag part
         const scalar EValsProd = ((EVals.x()*EVals.y()*EVals.z()).real());
-        cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant);
+        cmp("# Product of eigenvalues = det(T):", EValsProd, determinant, 1e-8);
     }
 
     {
@@ -929,7 +925,7 @@ void test_eigenvalues
         {
             EValsSum += val.real();
         }
-        cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace);
+        cmp("# Sum of eigenvalues = trace(T):", EValsSum, trace);
     }
 }
 
@@ -959,7 +955,7 @@ void test_characteristic_equation
 
         for (const auto x : X)
         {
-            cmp("  (sT & EVec - EVal*EVec) = 0:", mag(x), 0.0, 1e-8, 1e-6);
+            cmp("  (sT & EVec - EVal*EVec) = 0:", mag(x), 0.0, 1e-5);
         }
     }
 }
@@ -1056,14 +1052,13 @@ int main()
     }
 
     {
-        Info<< nl << "    ## Test eigen functions by a zero tensor: ##"<< nl;
+        Info<< nl << "    ## A zero tensor: ##"<< nl;
         const tensor zeroT(Zero);
         test_eigen_funcs(zeroT);
     }
     {
         Info<< nl
-            << "    ## Test eigen functions by a skew-symmetric tensor"
-            << " consisting of no-real eigenvalues: ##"
+            << "    ## A skew-symmetric tensor with no-real eigenvalues: ##"
             << nl;
         const tensor T
         (
@@ -1075,7 +1070,7 @@ int main()
     }
     {
         Info<< nl
-            << "    ## Test eigen functions by a stiff tensor: ##"
+            << "    ## A stiff tensor: ##"
             << nl;
         const tensor stiff
         (
@@ -1090,6 +1085,89 @@ int main()
         const bool testProd = false;
         test_eigen_funcs(stiff, testProd);
     }
+    {
+        Info<< nl
+            << "    ## Random tensor with tiny off-diag elements: ##"
+            << nl;
+
+        const List<scalar> epsilons
+        ({
+            0, SMALL, Foam::sqrt(SMALL), sqr(SMALL), Foam::cbrt(SMALL),
+            -SMALL, -Foam::sqrt(SMALL), -sqr(SMALL), -Foam::cbrt(SMALL)
+        });
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            for (const auto& eps : epsilons)
+            {
+                {
+                    tensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yz() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yx() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yx() = eps*rndGen.GaussNormal<scalar>();
+                    T.zx() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yx() = eps*rndGen.GaussNormal<scalar>();
+                    T.zx() = eps*rndGen.GaussNormal<scalar>();
+                    T.zy() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor T(makeRandomContainer(rndGen));
+                    T.xy() = 0;
+                    T.xz() = eps*rndGen.GaussNormal<scalar>();
+                    T.yz() = 0;
+                    T.yx() = eps*rndGen.GaussNormal<scalar>();
+                    T.zx() = eps*rndGen.GaussNormal<scalar>();
+                    T.zy() = 0;
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor T(makeRandomContainer(rndGen));
+                    T.xy() = eps;
+                    T.xz() = eps;
+                    T.yz() = eps;
+                    T.yx() = eps;
+                    T.zx() = eps;
+                    T.zy() = eps;
+                    test_eigen_funcs(T);
+                }
+            }
+        }
+    }
 
 
     if (nFail_)
diff --git a/applications/test/Tensor2D/Test-Tensor2D.C b/applications/test/Tensor2D/Test-Tensor2D.C
index b504e502fe..6563a8cb9c 100644
--- a/applications/test/Tensor2D/Test-Tensor2D.C
+++ b/applications/test/Tensor2D/Test-Tensor2D.C
@@ -86,11 +86,11 @@ typename std::enable_if
     const word& msg,
     const Type& x,
     const Type& y,
-    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
-    const scalar absTol = 0      //<! useful for cmps near zero
+    const scalar absTol = 0,     //<! useful for cmps near zero
+    const scalar relTol = 1e-8   //<! are values the same within 8 decimals
 )
 {
-    Info<< msg << x << endl;
+    Info<< msg << x << "?=" << y << endl;
 
     unsigned nFail = 0;
 
@@ -124,11 +124,11 @@ typename std::enable_if
     const word& msg,
     const Type& x,
     const Type& y,
-    const scalar relTol = 1e-8,
-    const scalar absTol = 0
+    const scalar absTol = 0,
+    const scalar relTol = 1e-8
 )
 {
-    Info<< msg << x << endl;
+    Info<< msg << x << "?=" << y << endl;
 
     unsigned nFail = 0;
 
@@ -725,7 +725,7 @@ void test_eigenvalues(const tensor2D& T, const Vector2D<complex>& EVals)
         // In case of complex EVals, the production is effectively scalar
         // due to the (complex*complex conjugate) results in zero imag part
         const scalar EValsProd = ((EVals.x()*EVals.y()).real());
-        cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant);
+        cmp("# Product of eigenvalues = det(T):", EValsProd, determinant, 1e-7);
     }
 
     {
@@ -737,7 +737,7 @@ void test_eigenvalues(const tensor2D& T, const Vector2D<complex>& EVals)
         {
             EValsSum += val.real();
         }
-        cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace, 1e-8, 1e-8);
+        cmp("# Sum of eigenvalues = trace(T):", EValsSum, trace, 1e-8);
     }
 }
 
@@ -767,7 +767,7 @@ void test_characteristic_equation
 
         for (const auto x : X)
         {
-            cmp("  (Tc & EVec - EVal*EVec) = 0:", mag(x), 0.0, 1e-8, 1e-6);
+            cmp("  (Tc & EVec - EVal*EVec) = 0:", mag(x), 0.0, 1e-5);
         }
     }
 }
@@ -861,14 +861,13 @@ int main(int argc, char *argv[])
     }
 
     {
-        Info<< nl << "    ## Test eigen functions by a zero tensor2D: ##"<< nl;
+        Info<< nl << "    ## A zero tensor2D: ##"<< nl;
         const tensor2D zeroT(Zero);
         test_eigen_funcs(zeroT);
     }
     {
         Info<< nl
-            << "    ## Test eigen functions by a tensor2D consisting of"
-            << " repeated eigenvalues: ##"
+            << "    ## A tensor2D with repeated eigenvalues: ##"
             << nl;
         const tensor2D T
         (
@@ -879,8 +878,7 @@ int main(int argc, char *argv[])
     }
     {
         Info<< nl
-            << "    ## Test eigen functions by a skew-symmetric tensor2D"
-            << " consisting of no-real eigenvalues: ##"
+            << "    ## A skew-symmetric tensor2D with no-real eigenvalues: ##"
             << nl;
         const tensor2D T
         (
@@ -891,7 +889,7 @@ int main(int argc, char *argv[])
     }
     {
         Info<< nl
-            << "    ## Test eigen functions by a stiff tensor2D: ##"
+            << "    ## A stiff tensor2D: ##"
             << nl;
         const tensor2D stiff
         (
@@ -900,6 +898,68 @@ int main(int argc, char *argv[])
         );
         test_eigen_funcs(stiff);
     }
+    {
+        Info<< nl
+            << "    ## Random tensor2D with tiny off-diag elements: ##"
+            << nl;
+
+        const List<scalar> epsilons
+        ({
+            0, SMALL, Foam::sqrt(SMALL), sqr(SMALL), Foam::cbrt(SMALL),
+            -SMALL, -Foam::sqrt(SMALL), -sqr(SMALL), -Foam::cbrt(SMALL)
+        });
+
+        for (label i = 0; i < numberOfTests; ++i)
+        {
+            for (const auto& eps : epsilons)
+            {
+                {
+                    tensor2D T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor2D T(makeRandomContainer(rndGen));
+                    T.yx() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor2D T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    T.yx() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor2D T(makeRandomContainer(rndGen));
+                    T.xy() = 0;
+                    T.yx() = eps*rndGen.GaussNormal<scalar>();
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor2D T(makeRandomContainer(rndGen));
+                    T.xy() = eps*rndGen.GaussNormal<scalar>();
+                    T.yx() = 0;
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor2D T(makeRandomContainer(rndGen));
+                    T.xy() = eps;
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor2D T(makeRandomContainer(rndGen));
+                    T.yx() = eps;
+                    test_eigen_funcs(T);
+                }
+                {
+                    tensor2D T(makeRandomContainer(rndGen));
+                    T.xy() = eps;
+                    T.yx() = eps;
+                    test_eigen_funcs(T);
+                }
+            }
+        }
+    }
 
 
     {
@@ -963,6 +1023,7 @@ int main(int argc, char *argv[])
         }
     }
 
+
     if (nFail_)
     {
         Info<< nl << "        #### "
diff --git a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
index b7a33e05e8..e8212aed3c 100644
--- a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
+++ b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
@@ -93,18 +93,27 @@ const Foam::symmTensor Foam::symmTensor::I
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Foam::vector Foam::eigenValues(const symmTensor& t)
+Foam::vector Foam::eigenValues(const symmTensor& T)
 {
+    // Return ascending diagonal if T is effectively diagonal tensor
+    if ((sqr(T.xy()) + sqr(T.xz()) + sqr(T.yz())) < ROOTSMALL)
+    {
+        vector eVals(T.diag());
+
+        std::sort(eVals.begin(), eVals.end());
+
+        return eVals;
+    }
+
     // Coefficients of the characteristic cubic polynomial (a = 1)
-    const scalar b =
-      - t.xx() - t.yy() - t.zz();
+    const scalar b = - T.xx() - T.yy() - T.zz();
     const scalar c =
-        t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
-      - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz();
+        T.xx()*T.yy() + T.xx()*T.zz() + T.yy()*T.zz()
+      - T.xy()*T.yx() - T.yz()*T.zy() - T.zx()*T.xz();
     const scalar d =
-      - t.xx()*t.yy()*t.zz()
-      - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx()
-      + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx();
+      - T.xx()*T.yy()*T.zz()
+      - T.xy()*T.yz()*T.zx() - T.xz()*T.zy()*T.yx()
+      + T.xx()*T.yz()*T.zy() + T.yy()*T.zx()*T.xz() + T.zz()*T.xy()*T.yx();
 
     // Determine the roots of the characteristic cubic polynomial
     Roots<3> roots(cubicEqn(1, b, c, d).roots());
@@ -123,26 +132,17 @@ Foam::vector Foam::eigenValues(const symmTensor& t)
             case roots::posInf:
             case roots::negInf:
             case roots::nan:
-                FatalErrorInFunction
-                    << "Eigenvalue computation fails for symmTensor: " << t
+                WarningInFunction
+                    << "Eigenvalue computation fails for symmTensor: " << T
                     << "due to the non-real root = " << roots[i]
-                    << exit(FatalError);
+                    << endl;
+                eVals[i] = roots[i];
+                break;
         }
     }
 
     // Sort the eigenvalues into ascending order
-    if (eVals.x() > eVals.y())
-    {
-        Swap(eVals.x(), eVals.y());
-    }
-    if (eVals.y() > eVals.z())
-    {
-        Swap(eVals.y(), eVals.z());
-    }
-    if (eVals.x() > eVals.y())
-    {
-        Swap(eVals.x(), eVals.y());
-    }
+    std::sort(eVals.begin(), eVals.end());
 
     return eVals;
 }
@@ -159,87 +159,89 @@ Foam::vector Foam::eigenVector
     // Construct the characteristic equation system for this eigenvalue
     const tensor A(T - eVal*I);
 
-    // Determinants of the 2x2 sub-matrices used to find the eigenvectors
-    // Sub-determinants for a unique eigenvenvalue
-    scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy();
-    scalar sd1 = A.zz()*A.xx() - A.zx()*A.xz();
-    scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx();
-    scalar magSd0 = mag(sd0);
-    scalar magSd1 = mag(sd1);
-    scalar magSd2 = mag(sd2);
-
-    // Evaluate the eigenvector using the largest sub-determinant
-    if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
     {
-        const vector eVec
-        (
-            1,
-            (A.yz()*A.zx() - A.zz()*A.yx())/sd0,
-            (A.zy()*A.yx() - A.yy()*A.zx())/sd0
-        );
+        // Determinants of the 2x2 sub-matrices used to find the eigenvectors
+        // Sub-determinants for a unique eigenvenvalue
+        const scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy();
+        const scalar sd1 = A.zz()*A.xx() - A.zx()*A.xz();
+        const scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx();
+        const scalar magSd0 = mag(sd0);
+        const scalar magSd1 = mag(sd1);
+        const scalar magSd2 = mag(sd2);
 
-        #ifdef FULLDEBUG
-        if (mag(eVec) < SMALL)
+        // Evaluate the eigenvector using the largest sub-determinant
+        if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
         {
-            FatalErrorInFunction
-                << "Eigenvector magnitude should be non-zero:"
-                << "mag(eigenvector) = " << mag(eVec)
-                << abort(FatalError);
+            const vector eVec
+            (
+                1,
+                (A.yz()*A.zx() - A.zz()*A.yx())/sd0,
+                (A.zy()*A.yx() - A.yy()*A.zx())/sd0
+            );
+
+            #ifdef FULLDEBUG
+            if (mag(eVec) < SMALL)
+            {
+                FatalErrorInFunction
+                    << "Eigenvector magnitude should be non-zero:"
+                    << "mag(eigenvector) = " << mag(eVec)
+                    << abort(FatalError);
+            }
+            #endif
+
+            return eVec/mag(eVec);
         }
-        #endif
-
-        return eVec/mag(eVec);
-    }
-    else if (magSd1 >= magSd2 && magSd1 > SMALL)
-    {
-        const vector eVec
-        (
-            (A.xz()*A.zy() - A.zz()*A.xy())/sd1,
-            1,
-            (A.zx()*A.xy() - A.xx()*A.zy())/sd1
-        );
-
-        #ifdef FULLDEBUG
-        if (mag(eVec) < SMALL)
+        else if (magSd1 >= magSd2 && magSd1 > SMALL)
         {
-            FatalErrorInFunction
-                << "Eigenvector magnitude should be non-zero:"
-                << "mag(eigenvector) = " << mag(eVec)
-                << abort(FatalError);
+            const vector eVec
+            (
+                (A.xz()*A.zy() - A.zz()*A.xy())/sd1,
+                1,
+                (A.zx()*A.xy() - A.xx()*A.zy())/sd1
+            );
+
+            #ifdef FULLDEBUG
+            if (mag(eVec) < SMALL)
+            {
+                FatalErrorInFunction
+                    << "Eigenvector magnitude should be non-zero:"
+                    << "mag(eigenvector) = " << mag(eVec)
+                    << abort(FatalError);
+            }
+            #endif
+
+            return eVec/mag(eVec);
         }
-        #endif
-
-        return eVec/mag(eVec);
-    }
-    else if (magSd2 > SMALL)
-    {
-        const vector eVec
-        (
-            (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
-            (A.yx()*A.xz() - A.xx()*A.yz())/sd2,
-            1
-        );
-
-        #ifdef FULLDEBUG
-        if (mag(eVec) < SMALL)
+        else if (magSd2 > SMALL)
         {
-            FatalErrorInFunction
-                << "Eigenvector magnitude should be non-zero:"
-                << "mag(eigenvector) = " << mag(eVec)
-                << abort(FatalError);
-        }
-        #endif
+            const vector eVec
+            (
+                (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
+                (A.yx()*A.xz() - A.xx()*A.yz())/sd2,
+                1
+            );
 
-        return eVec/mag(eVec);
+            #ifdef FULLDEBUG
+            if (mag(eVec) < SMALL)
+            {
+                FatalErrorInFunction
+                    << "Eigenvector magnitude should be non-zero:"
+                    << "mag(eigenvector) = " << mag(eVec)
+                    << abort(FatalError);
+            }
+            #endif
+
+            return eVec/mag(eVec);
+        }
     }
 
     // Sub-determinants for a repeated eigenvalue
-    sd0 = A.yy()*standardBasis1.z() - A.yz()*standardBasis1.y();
-    sd1 = A.zz()*standardBasis1.x() - A.zx()*standardBasis1.z();
-    sd2 = A.xx()*standardBasis1.y() - A.xy()*standardBasis1.x();
-    magSd0 = mag(sd0);
-    magSd1 = mag(sd1);
-    magSd2 = mag(sd2);
+    const scalar sd0 = A.yy()*standardBasis1.z() - A.yz()*standardBasis1.y();
+    const scalar sd1 = A.zz()*standardBasis1.x() - A.zx()*standardBasis1.z();
+    const scalar sd2 = A.xx()*standardBasis1.y() - A.xy()*standardBasis1.x();
+    const scalar magSd0 = mag(sd0);
+    const scalar magSd1 = mag(sd1);
+    const scalar magSd2 = mag(sd2);
 
     // Evaluate the eigenvector using the largest sub-determinant
     if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
diff --git a/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.C b/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.C
index 7ff6ade25e..5eec5b82d0 100644
--- a/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.C
+++ b/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.C
@@ -88,12 +88,18 @@ const Foam::symmTensor2D Foam::symmTensor2D::I
 
 Foam::vector2D Foam::eigenValues(const symmTensor2D& T)
 {
+    // Return diagonal if T is effectively diagonal tensor
+    if (sqr(T.xy()) < ROOTSMALL)
+    {
+        return vector2D(T.diag());
+    }
+
     //(K:Eqs. 3.2-3.3)
     const scalar skewTrace = T.xx() - T.yy();
     const scalar trace = tr(T);
     const scalar gap = sign(skewTrace)*hypot(skewTrace, 2.0*T.xy());
 
-    return vector2D (0.5*(trace + gap), 0.5*(trace - gap));
+    return vector2D(0.5*(trace + gap), 0.5*(trace - gap));
 }
 
 
diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
index 69471de43f..776c11083c 100644
--- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
+++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
@@ -74,22 +74,35 @@ const Foam::tensor Foam::tensor::I
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Foam::Vector<Foam::complex> Foam::eigenValues(const tensor& t)
+Foam::Vector<Foam::complex> Foam::eigenValues(const tensor& T)
 {
+    // Return diagonal if T is effectively diagonal tensor
+    if
+    (
+        (
+            sqr(T.xy()) + sqr(T.xz()) + sqr(T.yz())
+          + sqr(T.yx()) + sqr(T.zx()) + sqr(T.zy())
+        ) < ROOTSMALL
+    )
+    {
+        return Vector<complex>
+        (
+            complex(T.xx()), complex(T.yy()), complex(T.zz())
+        );
+    }
+
     // Coefficients of the characteristic cubic polynomial (a = 1)
-    const scalar b =
-      - t.xx() - t.yy() - t.zz();
+    const scalar b = - T.xx() - T.yy() - T.zz();
     const scalar c =
-        t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
-      - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz();
+        T.xx()*T.yy() + T.xx()*T.zz() + T.yy()*T.zz()
+      - T.xy()*T.yx() - T.yz()*T.zy() - T.zx()*T.xz();
     const scalar d =
-      - t.xx()*t.yy()*t.zz()
-      - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx()
-      + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx();
+      - T.xx()*T.yy()*T.zz()
+      - T.xy()*T.yz()*T.zx() - T.xz()*T.zy()*T.yx()
+      + T.xx()*T.yz()*T.zy() + T.yy()*T.zx()*T.xz() + T.zz()*T.xy()*T.yx();
 
     // Determine the roots of the characteristic cubic polynomial
     const Roots<3> roots(cubicEqn(1, b, c, d).roots());
-
     // Check the root types
     bool isComplex = false;
     forAll(roots, i)
@@ -102,10 +115,10 @@ Foam::Vector<Foam::complex> Foam::eigenValues(const tensor& t)
             case roots::posInf:
             case roots::negInf:
             case roots::nan:
-                FatalErrorInFunction
-                    << "Eigenvalue computation fails for tensor: " << t
+                WarningInFunction
+                    << "Eigenvalue computation fails for tensor: " << T
                     << "due to the not-a-number root = " << roots[i]
-                    << exit(FatalError);
+                    << endl;
             case roots::real:
                 break;
         }
diff --git a/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C b/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C
index 4cfca5d440..93911913f1 100644
--- a/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C
+++ b/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C
@@ -94,6 +94,12 @@ Foam::Vector2D<Foam::complex> Foam::eigenValues(const tensor2D& T)
     const scalar c = T.yx();
     const scalar d = T.yy();
 
+    // Return diagonal if T is effectively diagonal tensor
+    if ((sqr(b) + sqr(c)) < ROOTSMALL)
+    {
+        return Vector2D<complex>(complex(a), complex(d));
+    }
+
     const scalar trace = a + d;
 
     // (JLM:p. 2246)
@@ -126,7 +132,7 @@ Foam::Vector2D<Foam::complex> Foam::eigenValues(const tensor2D& T)
         );
 
         // Sort the eigenvalues into ascending order
-        if (mag(eVals.x()) > mag(eVals.y()))
+        if (eVals.x().real() > eVals.y().real())
         {
             Swap(eVals.x(), eVals.y());
         }
@@ -154,9 +160,48 @@ Foam::Vector2D<Foam::complex> Foam::eigenVector
     const Vector2D<complex>& standardBasis
 )
 {
-    // (K:p. 47-48)
+    // Construct the linear system for this eigenvalue
+    const Tensor2D<complex> A
+    (
+        complex(T.xx()) - eVal,  complex(T.xy()),
+        complex(T.yx()),         complex(T.yy()) - eVal
+    );
+
     // Evaluate the eigenvector using the largest divisor
-    if (mag(T.yx()) > mag(T.xy()) && mag(T.yx()) > VSMALL)
+    if (mag(A.yy()) > mag(A.xx()) && mag(A.yy()) > SMALL)
+    {
+        Vector2D<complex> eVec(complex(1), -A.yx()/A.yy());
+
+        #ifdef FULLDEBUG
+        if (mag(eVec) < SMALL)
+        {
+            FatalErrorInFunction
+                << "Eigenvector magnitude should be non-zero:"
+                << "mag(eigenvector) = " << mag(eVec)
+                << abort(FatalError);
+        }
+        #endif
+
+        return eVec/mag(eVec);
+    }
+    else if (mag(A.xx()) > SMALL)
+    {
+        Vector2D<complex> eVec(-A.xy()/A.xx(), complex(1));
+
+        #ifdef FULLDEBUG
+        if (mag(eVec) < SMALL)
+        {
+            FatalErrorInFunction
+                << "Eigenvector magnitude should be non-zero:"
+                << "mag(eigenvector) = " << mag(eVec)
+                << abort(FatalError);
+        }
+        #endif
+
+        return eVec/mag(eVec);
+    }
+    // (K:p. 47-48)
+    else if (mag(T.yx()) > mag(T.xy()) && mag(T.yx()) > SMALL)
     {
         const Vector2D<complex> eVec(eVal - T.yy(), complex(T.yx()));
 
@@ -172,7 +217,7 @@ Foam::Vector2D<Foam::complex> Foam::eigenVector
 
         return eVec/mag(eVec);
     }
-    else if (mag(T.xy()) > VSMALL)
+    else if (mag(T.xy()) > SMALL)
     {
         const Vector2D<complex> eVec(complex(T.xy()), eVal - T.xx());
 
@@ -189,7 +234,8 @@ Foam::Vector2D<Foam::complex> Foam::eigenVector
         return eVec/mag(eVec);
     }
 
-    return standardBasis;
+    // Repeated eigenvalue
+    return Vector2D<complex>(-standardBasis.y(), standardBasis.x());
 }
 
 
diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C
index fc69445e97..17732e91dc 100644
--- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C
+++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C
@@ -39,6 +39,16 @@ Foam::Roots<3> Foam::cubicEqn::roots() const
     const scalar c = this->c();
     const scalar d = this->d();
 
+    #ifdef FULLDEBUG
+    Info<< "#DEBUG#" << nl
+        << "Coefficients of the characteristic cubic polynomial:" << nl
+        << "a = " << a << nl
+        << "b = " << b << nl
+        << "c = " << c << nl
+        << "d = " << d << nl
+        << "#######" << endl;
+    #endif
+
     // Check the leading term in the cubic eqn exists
     if (mag(a) < VSMALL)
     {
@@ -50,15 +60,27 @@ Foam::Roots<3> Foam::cubicEqn::roots() const
     const scalar p = -(fma(-a, c, w) + fma(b, b/3.0, -w));
     const scalar q = b*b*b*scalar(2)/27 - b*c*a/3 + d*a*a;
     const scalar numDiscr = p*p*p/27 + q*q/4;
-    const scalar discr = (mag(numDiscr) > SMALL) ? numDiscr : 0;
+    const scalar discr = (mag(numDiscr) > sqr(SMALL)) ? numDiscr : 0;
 
     // Determine the number and types of the roots
     const bool threeReal = discr < 0;
     const bool oneRealTwoComplex = discr > 0;
     const bool twoReal = //p != 0; && discr == 0;
-        (mag(p) > VSMALL) && !(threeReal || oneRealTwoComplex);
+        (mag(p) > sqrt(SMALL)) && !(threeReal || oneRealTwoComplex);
     // const bool oneReal = p == 0 && discr == 0;
 
+    #ifdef FULLDEBUG
+    Info<< "#DEBUG#" << nl
+        << "Numerical discriminant:" << tab << numDiscr << nl
+        << "Adjusted discriminant:" << tab << discr << nl
+        << "Number and types of the roots:" << nl
+        << "threeReal = " << threeReal << nl
+        << "oneRealTwoComplex = " << oneRealTwoComplex << nl
+        << "twoReal = " << twoReal << nl
+        << "oneReal = " << !(threeReal || oneRealTwoComplex || twoReal) << nl
+        << "#######" << endl;
+    #endif
+
     static const scalar sqrt3 = sqrt(3.0);
 
     scalar x = 0;
@@ -127,6 +149,12 @@ Foam::Roots<3> Foam::cubicEqn::roots() const
         return Roots<3>(r.type(0), r[0]);
     }
 
+    #if FULLDEBUG
+    Info<< "#DEBUG#" << nl
+        << "x = " << x << nl
+        << "#######" << endl;
+    #endif
+
     return
         Roots<3>
         (