Timo,
Merge request done. Hope it will make it to DuMux 3.4.
Best regards,
Dmitry
On 30.04.2021 00:10, Dmitry Pavlov wrote:
Timo,
I will do this next week.
Best regards,
Dmitry
On 29.04.2021 17:22, Timo Koch wrote:
Hi Dmitry,
thanks for looking into this. It would be very nice if you could
provide an MR with the proposed fix
together with some unit tests that test corner cases (in
test/common/test_math.cc).
Timo
On 29. Apr 2021, at 16:05, Dmitry Pavlov <dmitry.pav...@outlook.com>
wrote:
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
DuMux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________
DuMux mailing list
DuMux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________
DuMux mailing list
DuMux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
