Fix of a corner case in inverseCubicPolynomial + test for that corner case

dumux/common/math.hh

@@ -330,8 +330,14 @@ int invertCubicPolynomial(SolContainer *sol,
         // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
         using std::cbrt;
         using std::sqrt;
-        Scalar u = cbrt(-q/2 + sqrt(wDisc));
 
+        // 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. 
+        const Scalar u = [&]{
+            return q > 0 ? cbrt(-0.5*q - sqrt(wDisc)) : cbrt(-0.5*q + sqrt(wDisc));
+        }();
         // at this point, u != 0 since p^3 = 0 is necessary in order
         // for u = 0 to hold, so
         sol[0] = u - p/(3*u) - b/3;

test/common/math/test_math.cc

@@ -235,4 +235,47 @@ int main()
             DUNE_THROW(Dune::Exception, "[linearRegreesion] Wrong slope " << slope << ", should be 2.0.");
     }
 
+    //////////////////////////////////////////////////////////////////
+    ///// Dumux::invertCubicPolynomial ///////////////////////////////
+    //////////////////////////////////////////////////////////////////
+
+    std::array<double, 3> roots{};
+
+    // One corner case
+    int nroots = Dumux::invertCubicPolynomial(roots.data(),
+                                              1.0,
+                                              -0.96165703943410097,
+                                              0.30826068470787077,
+                                              0.01340255221155587);
+    if (nroots != 1)
+        DUNE_THROW(Dune::Exception, "Expecting one root");
+
+    using std::isfinite;
+    
+    if (!isfinite(roots[0]))
+        DUNE_THROW(Dune::Exception, "Expecting finite root");
+    if (!Dune::FloatCmp::eq(roots[0], -3.863448718244389e-02, 1e-14))
+        DUNE_THROW(Dune::Exception, "Root of cubic equation does not match reference");
+
+
+    // One simple case
+    std::array<double, 3> known_roots{-10.1, 0.001, 1000.77};
+
+    nroots = Dumux::invertCubicPolynomial(
+        roots.data(),
+        1.0,
+        -known_roots[0] - known_roots[1] - known_roots[2],
+        known_roots[0] * known_roots[1] + known_roots[0] * known_roots[2] + known_roots[1] * known_roots[2],
+        -known_roots[0] * known_roots[1] * known_roots[2]);
+
+    if (nroots != 3)
+        DUNE_THROW(Dune::Exception, "Expecting three roots");
+
+    for (int i = 0; i < 3; i++)
+    {
+        if (!isfinite(roots[i]))
+            DUNE_THROW(Dune::Exception, "Expecting finite root");
+        if (!Dune::FloatCmp::eq(roots[i], known_roots[i], 1e-10))
+            DUNE_THROW(Dune::Exception, "Root of cubic equation does not match reference");
+    }
 }