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 <[email protected]> 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 > [email protected] > https://listserv.uni-stuttgart.de/mailman/listinfo/dumux _______________________________________________ DuMux mailing list [email protected] https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
