Re: [DuMux] Numerical bug in invertCubicPolynomial

[Timo Koch]

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  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] Numerical bug in invertCubicPolynomial

[Dmitry Pavlov]

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] [DuMuX] DuMuX 3.4 Release Schedule and Important Information

[Wu, Hanchuan]

Dear DuMuX community,

Dumux 3.4 is scheduled to be released on May 18th 2021!

Please consider the following dates while developing new features and adding 
merge requests to Dumux:

== Soft Feature Freeze: 28. April 2021 (3 weeks prior to the release) ==
- All major changes that should be included and are not yet committed should be 
announced. Add a comment to 
Issue#1009
 describing the feature and the areas affected by the change.
- Major changes not announced before the soft freeze will be reserved for the 
next release.
- The representative assigned to each major task will report on their remaining 
tasks.

== Hard Feature Freeze: 5th. May 2021 (2 weeks prior to the release) ==
- All major changes must be committed before the hard feature freeze.
- Large invasive fixes are only accepted after consultation with the release 
manager (small fixes are ok).
- Try to get your bugs fixed, your tests included, and your documentation 
written before the hard feature freeze.
- The representative assigned to each major task will give a final report on 
what was finished and what must be postponed until the next release.

== Final Testing:  12. May 2021 (1 week prior to the release) ==
- No further commits are allowed, barring those performed by the release 
manager, or under the explicit authorisation of the release manager.
- Release candidate packages will be provided.
- Tests will be performed on various systems, compilers, and under different 
conditions. Everyone should at least test their own tests!

== DuMuX 3.4 Release: 19. May 2021 ==
- Celebrate! Virtual cold beverages and schnapps are available.

Looking forward to a great release, Hanchuan (Release manager DuMuX 3.4)

***
Hanchuan Wu
Universitaet Stuttgart
Lehrstuhl fuer Hydromechanik & Hydrosystemmodellierung
Pfaffenwaldring 61
70569 Stuttgart
Tel.: +49-711-685 69163
hanchuan...@iws.uni-stuttgart.de
https://www.iws.uni-stuttgart.de/lh2/
***
___
DuMux mailing list
DuMux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux

[DuMux] How to install dune-alugrid

[Deming Wang]

Hi experts,

I tried to install dune-alugrid in my dumux folder(fresh installation),
and executed the following command:
 ./dumux/bin/installexternal.sh all
then run duneproject:
[user@server1 dumux]$ ./dune-common/bin/duneproject
and end with errors:

WARNING: could not find module 'dune-uggrid',
   module is also unknown to pkg-config.
   Maybe you need to adjust PKG_CONFIG_PATH!
   'dune-uggrid' is suggested by dune-grid
Skipping 'dune-uggrid'!
ERROR: version mismatch.
   dune-alugrid requires 'dune-grid >= 2.8',
   but only 'dune-grid' = '2.7.1' is available.

I try to upgrade to dune-grid 2.8, but there is no dune-grid release yet.


Any instructions about how to install dune-alugrid?  Thanks.


Regards,

Deming

-- 
Deming Wang
E-mail:wangdem...@gmail.com
___
DuMux mailing list
DuMux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux