Hello,
I would like to propose a fix to invertCubicPolynomial
(dumux/common/math.hh) in the following block:
Scalar wDisc = q*q/4 + p*p*p/27;
if (wDisc >= 0) { // the positive discriminant case:
// 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));
// 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;
// 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);
return 1;
}
In fact, u can still be zero with a nonzero p due to floating-point
catastrophic cancellation.
Example where that happens:
a = 1
b = -0.96165703943410097
c = 0.30826068470787077
d = 0.01340255221155587
To avoid the division by zero (and in general to avoid the
ill-conditioned way of calculation), I would suggest to choose the other
root of the quadratic equation. Generally speaking, we want the root for
which the catastrophic cancellation will not occur.
if (q > 0)
u = cbrt(-q/2 - sqrt(wDisc));
else
u = cbrt(-q/2 + sqrt(wDisc));
We are happy with either root, so it should work. (If we need the other
root, too, for some reason, we can obtain it from the first one using
the Vieta's formula).
Best regards,
Dmitry
_______________________________________________
DuMux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
