Merge branch 'fix/invert-poly-doc-and-test' into 'master'

[math][invertCubicPoly] Improve doc. Add possibility to specify...

Closes #1025

See merge request !2599

dumux/common/math.hh

@@ -198,25 +198,30 @@ int invertQuadraticPolynomial(SolContainer &sol,
 template <class Scalar, class SolContainer>
 void invertCubicPolynomialPostProcess_(SolContainer &sol,
                                        int numSol,
-                                       Scalar a,
-                                       Scalar b,
-                                       Scalar c,
-                                       Scalar d)
+                                       Scalar a, Scalar b, Scalar c, Scalar d,
+                                       std::size_t numIterations = 1)
 {
-    // do one Newton iteration on the analytic solution if the
-    // precision is increased
-    for (int i = 0; i < numSol; ++i) {
-        Scalar x = sol[i];
-        Scalar fOld = d + x*(c + x*(b + x*a));
+    const auto eval = [&](auto x){ return d + x*(c + x*(b + x*a)); };
+    const auto evalDeriv = [&](auto x){ return c + x*(2*b + x*3*a); };
 
-        Scalar fPrime = c + x*(2*b + x*3*a);
-        if (fPrime == 0.0)
-            continue;
-        x -= fOld/fPrime;
+    // do numIterations Newton iterations on the analytic solution
+    // and update result if the precision is increased
+    for (int i = 0; i < numSol; ++i)
+    {
+        Scalar x = sol[i];
+        Scalar fCurrent = eval(x);
+        const Scalar fOld = fCurrent;
+        for (int j = 0; j < numIterations; ++j)
+        {
+            const Scalar fPrime = evalDeriv(x);
+            if (fPrime == 0.0)
+                break;
+            x -= fCurrent/fPrime;
+            fCurrent = eval(x);
+        }
 
-        Scalar fNew = d + x*(c + x*(b + x*a));
         using std::abs;
-        if (abs(fNew) < abs(fOld))
+        if (abs(fCurrent) < abs(fOld))
             sol[i] = x;
     }
 }
@@ -229,22 +234,24 @@ void invertCubicPolynomialPostProcess_(SolContainer &sol,
  * The polynomial is defined as
  * \f[ p(x) = a\; x^3 + + b\;x^3 + c\;x + d \f]
  *
- * This method teturns the number of solutions which are in the real
+ * This method returns the number of solutions which are in the real
  * numbers. The "sol" argument contains the real roots of the cubic
  * polynomial in order with the smallest root first.
  *
+ * \note The closer the roots are to each other the less
+ *       precise the inversion becomes. Increase number of post-processing iterations for improved results.
+ *
  * \param sol Container into which the solutions are written
- * \param a The coefficiont for the cubic term
- * \param b The coefficiont for the quadratic term
- * \param c The coefficiont for the linear term
- * \param d The coefficiont for the constant term
+ * \param a The coefficient for the cubic term
+ * \param b The coefficient for the quadratic term
+ * \param c The coefficient for the linear term
+ * \param d The coefficient for the constant term
+ * \param numPostProcessIterations The number of iterations to increase precision of the analytical result
  */
 template <class Scalar, class SolContainer>
 int invertCubicPolynomial(SolContainer *sol,
-                          Scalar a,
-                          Scalar b,
-                          Scalar c,
-                          Scalar d)
+                          Scalar a, Scalar b, Scalar c, Scalar d,
+                          std::size_t numPostProcessIterations = 1)
 {
     // reduces to a quadratic polynomial
     if (a == 0)
@@ -334,7 +341,7 @@ int invertCubicPolynomial(SolContainer *sol,
         // Choose the root that is safe against the loss of precision
         // that can cause wDisc to be equal to q*q/4 despite p != 0.
         // We do not want u to be zero in that case. Mathematically,
-        // we are happy with either root. 
+        // we are happy with either root.
         const Scalar u = [&]{
             return q > 0 ? cbrt(-0.5*q - sqrt(wDisc)) : cbrt(-0.5*q + sqrt(wDisc));
         }();
@@ -344,7 +351,7 @@ int invertCubicPolynomial(SolContainer *sol,
         // does not produce a division by zero. the remaining two
         // roots of u are rotated by +- 2/3*pi in the complex plane
         // and thus not considered here
-        invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d);
+        invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d, numPostProcessIterations);
         return 1;
     }
     else { // the negative discriminant case:
@@ -404,7 +411,7 @@ int invertCubicPolynomial(SolContainer *sol,
 
         // post process the obtained solution to increase numerical
         // precision
-        invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d);
+        invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d, numPostProcessIterations);
 
         // sort the result
         using std::sort;

test/common/math/test_math.cc

@@ -44,6 +44,7 @@
 #include <dune/common/dynvector.hh>
 
 #include <dumux/common/math.hh>
+#include <dumux/io/format.hh>
 
 namespace Test {
 
@@ -264,31 +265,43 @@ int main()
     std::uniform_real_distribution<> dis(-10.0, 10.0);
     auto threeRandomNumbers = [&]() mutable {
         const auto x = dis(gen);
-        auto y = dis(gen); while (x == y) { y = dis(gen); }
-        auto z = dis(gen); while (x == z || y == z) { z = dis(gen); }
+        auto y = dis(gen); while (Dune::FloatCmp::eq(x, y, 1e-10)) { y = dis(gen); }
+        auto z = dis(gen); while (Dune::FloatCmp::eq(x, z, 1e-10) || Dune::FloatCmp::eq(y, z, 1e-10)) { z = dis(gen); }
         return std::make_tuple(x, y, z);
     };
 
     const auto checkRoots = [](const std::array<double, 3>& roots,
                                const std::array<double, 3>& ref,
-                               int n)
+                               const std::array<double, 4>& coeff,
+                               int n, bool throwing = true) -> bool
     {
+        const auto refToString = [&](){
+            std::string refStr{};
+            for (int i = 0; i < n; i++)
+                refStr += Dumux::Fmt::format("{:.14e},", ref[i]);
+            return refStr;
+        };
+
         for (int i = 0; i < n; i++)
         {
             using std::isfinite;
             if (!isfinite(roots[i]))
                 DUNE_THROW(Dune::Exception, "Expecting finite root");
-            if (std::none_of(ref.begin(), ref.end(), [&](auto r){ return Dune::FloatCmp::eq(roots[i], r, 1e-10); }))
-                DUNE_THROW(Dune::Exception, "[invertCubicPolynomial] Root "
-                    << i << " of " << n << ": " << roots[i]
-                    << " does not match reference");
+            if (throwing && std::none_of(ref.begin(), ref.end(), [&](auto r){ return Dune::FloatCmp::eq(roots[i], r, 1e-8); }))
+                DUNE_THROW(Dune::Exception, "[invertCubicPolynomial] Root " << Dumux::Fmt::format(
+                        "{} of {}: {:.14e}, does not match reference roots = [{}].\n", i+1, n, roots[i], refToString())
+                        << Dumux::Fmt::format("Polynomial: {:+.5e}x³ {:+.5e}x² {:+.5e}x {:+.5e}\n", coeff[0], coeff[1], coeff[2], coeff[3]));
+            else
+                return false;
         }
+
+        return true;
     };
 
     std::array<double, 3> roots{};
     for (double scaling : {1.0, -1.0, 5.0})
     {
-        for (int i = 0; i < 1000; ++i)
+        for (int i = 0; i < 100000; ++i)
         {
             const auto [x, y, z] = threeRandomNumbers();
 
@@ -298,7 +311,12 @@ int main()
                 int numRoots = Dumux::invertCubicPolynomial(roots.data(), a, b, c, d);
                 if (numRoots != 3)
                     DUNE_THROW(Dune::Exception, "Expecting three roots! Got " << numRoots);
-                checkRoots(roots, {x, y, z}, numRoots);
+                if (!checkRoots(roots, {x, y, z}, {a, b, c, d}, numRoots, false))
+                {
+                    // Try to increase Newton iterations for increased precision and try again
+                    int numRoots = Dumux::invertCubicPolynomial(roots.data(), a, b, c, d, 10);
+                    checkRoots(roots, {x, y, z}, {a, b, c, d}, numRoots);
+                }
             }
 
             // test for two unique roots
@@ -307,7 +325,12 @@ int main()
                 int numRoots = Dumux::invertCubicPolynomial(roots.data(), 0.0, b, c, d);
                 if (numRoots != 2)
                     DUNE_THROW(Dune::Exception, "Expecting two roots! Got " << numRoots);
-                checkRoots(roots, {x, z}, numRoots);
+                if (!checkRoots(roots, {x, z}, {0.0, b, c, d}, numRoots, false))
+                {
+                    // Try to increase Newton iterations for increased precision and try again
+                    int numRoots = Dumux::invertCubicPolynomial(roots.data(), 0.0, b, c, d, 10);
+                    checkRoots(roots, {x, z}, {0.0, b, c, d}, numRoots);
+                }
             }
 
             // test for one unique root
@@ -316,7 +339,12 @@ int main()
                 int numRoots = Dumux::invertCubicPolynomial(roots.data(), 0.0, 0.0, c, d);
                 if (numRoots != 1)
                     DUNE_THROW(Dune::Exception, "Expecting one root! Got " << numRoots);
-                checkRoots(roots, {x}, numRoots);
+                if (!checkRoots(roots, {x}, {0.0, 0.0, c, d}, numRoots, false))
+                {
+                    // Try to increase Newton iterations for increased precision and try again
+                    int numRoots = Dumux::invertCubicPolynomial(roots.data(), 0.0, 0.0, c, d, 10);
+                    checkRoots(roots, {x}, {0.0, 0.0, c, d}, numRoots);
+                }
             }
         }
     }