View Issue Details

IDProjectCategoryView StatusLast Update
0002972OpenFOAMBugpublic2018-06-26 21:02
Reporterpeksa Assigned Tohenry  
PrioritylowSeverityminorReproducibilityalways
Status closedResolutionsuspended 
PlatformGNU/LinuxOSOtherOS Version(please specify)
Fixed in Versiondev 
Summary0002972: Stability / monotonicity check issue in the Seulex ODE solver
DescriptionIn commit 1ef6ed3b9f34766af8c835c9e5c6850ba0f5c668 (OF-5.x) the stability / monotonicity check of the seulex algorithm (in seulex.C, function seul() ) was modified with an aim to prevent possible overflow situations.

However, currently, the monotonicity check (OF-5.x code lines 156-170) is not following the original implementation by Hairer et al. http://www.unige.ch/~hairer/prog/stiff/seulex.f code lines 1110-1115. In particular the denominator
"const scalar denom = min(1, dy1 + SMALL);"
is originally defined as max(1, dy1)" as it was also done in the earlier OpenFOAM versions.

When trying to reach long integration steps, the current implementation results overhead because typically dy1 is a large number and thus "denom" becomes equal to unity due to the min function, resulting in algorithm rather easily discarding that time step in the following if-closure:

if (mag(dy_[i]) > scale[i]*denom)
{
    theta_ = 1;
    return false;
}

However, this seems not to be a big problem in CFD applications where we anyway have a rather small initial time step and thus less discarded steps due to monotonicity check.

I'm not so sure is the original (Hairer et al.) approach for this correct either but I have my doubt of the generality of the current OF approach as well. Monotonicity check is a tricky thing with extrapolation methods, especially when trying to be eager with long ODE time steps.

Furthermore, my experience with seulex does not include any problems regarding overflow in this monotonicity check (old implementation). When I apply it to my chemistry computations I use a third-party reaction rate and jacobian computations and have never had any issues with any mechanism, so the found overflows most probably rise due to a not well-defined ODE problem found in certain cells of multiphysics cases.


Steps To ReproduceI will provide an example of the issue later if needed.
TagsODE, seulex

Activities

henry

2018-06-05 11:29

manager   ~0009704

It is not clear from your description what change you are proposing, could you clarify?

peksa

2018-06-06 07:56

reporter   ~0009705

I'd revert this stability check part back to the old version (before the commit mentioned above) which is consistent with the original Hairer's implementation. Particularly, the min function should be changed to max function in the examples above.

henry

2018-06-06 08:50

manager   ~0009706

If I do that the cases which failed before the change will fail again. I made the change specifically to resolve overflow issues when solving very stiff problems. However, I can check the min/max change on a range of cases.

henry

2018-06-06 11:41

manager   ~0009707

> is originally defined as max(1, dy1)" as it was also done in the earlier OpenFOAM versions.

Which earlier versions are you referring to? Could provide which version, file and line number?

peksa

2018-06-06 12:21

reporter   ~0009708

Ahh, damn, I somehow got mixed up with the OF original repository versions and my own old implementations for OF-2.4.x. I'm sorry for misleading you.

Indeed, the original old OF-2.3.x-2.4.x implementations have that min-function there as well and I had fixed it to max on my own when I was going through the Hairer's implementation and did not do a bug report at that time.

Furthermore, Authors of https://www.sciencedirect.com/science/article/pii/S001021801630267X#sec0024
have published an OF based seulex algorithm (see supplementary material or the attached pdf here) where they apply the following stability check in the seul() function. Otherwise their published piece of code shows just how to include 3rd-party sparse linear algebra into the seulex code.


for (label i=0; i<n_; i++)
{
    // check magnitude of change in solution vector, compared to VGREAT
    if (mag(dy_[i]) > sqrt(VGREAT) ){ dy2 = VGREAT; break;}
    dy2 += sqr(dy_[i]/scale[i]);
}
dy2 = sqrt(dy2);
// check convergence of Newton method: apply monotonicity check
theta_ = dy2/max(1.0, dy1 + SMALL);

The above is in line what I was talking previously and comparable to the Hairer's original piece of Fortran code. In some point I omitted the stability check from that publication (VGREAT comparison) and have not had any problems regarding overflow nor other crashes.

Once again, I'm sorry for the wrong information in the beginning.
seulex_Imren.pdf (777,261 bytes)

henry

2018-06-06 12:36

manager   ~0009709

Thanks for the clarification. Unfortunately the authors of the paper you cite did not report the min/max issues either otherwise this could have been resolved some time ago.

I am testing the change

            scalar dy2 = 0;
            for (label i=0; i<n_; i++)
            {
                dy2 += sqr(dy_[i]/scale[i]);
            }
            dy2 = sqrt(dy2);
            theta_ = dy2/max(1, dy1);

and in this form it is not clear that the overflow test is necessary, perhaps the authors of the paper can provide the reasoning for the overflow test they implemented.

henry

2018-06-06 12:49

manager   ~0009710

I can generate overflow for VERY stiff problems and reinstated my overflow test which I think is preferable to the one in the paper you cite:

            const scalar denom = max(1, dy1);
            scalar dy2 = 0;
            for (label i=0; i<n_; i++)
            {
                // Test of dy_[i] to avoid overflow
                if (mag(dy_[i]) > scale[i]*denom)
                {
                    theta_ = 1;
                    return false;
                }

                dy2 += sqr(dy_[i]/scale[i]);
            }
            dy2 = sqrt(dy2);
            theta_ = dy2/denom;

If this is also OK for your cases and you do not see any other issues with it I will push it into OpenFOAM-dev and 5.x

peksa

2018-06-07 07:49

reporter   ~0009715

I can confirm that your proposed implementation works for my cases as well.

henry

2018-06-07 08:04

manager   ~0009716

Resolved in OpenFOAM-5.x by commit 7f7d351b741bf6406366a043cac98de56d2d44dd
Resolved in OpenFOAM-dev by commit 873d2c9c23a617432f96ca5cb9ffe6c02661e9c1

henry

2018-06-11 11:16

manager   ~0009737

With this change the tutorials/lagrangian/sprayFoam/aachenBomb case no longer runs:

--> FOAM Warning :
    From function virtual void Foam::seulex::solve(Foam::scalar&, Foam::scalarField&, Foam::ODESolver::stepState&) const
    in file ODESolvers/seulex/seulex.C at line 289
    step size underflow :7.89336e-41
.
.
.

Do you know how to resolve this issue? If not I will have to revert the change.

peksa

2018-06-12 14:45

reporter   ~0009742

I will look into this tomorrow when I get back to my computer.

henry

2018-06-22 12:08

manager   ~0009811

Do you have any suggestions for alternative formulations which are reliable for all cases? If not should I keep the original formulation and close this report?

henry

2018-06-26 21:02

manager   ~0009833

Pending feedback from reporter

Issue History

Date Modified Username Field Change
2018-06-05 11:07 peksa New Issue
2018-06-05 11:07 peksa Tag Attached: ODE
2018-06-05 11:07 peksa Tag Attached: seulex
2018-06-05 11:29 henry Note Added: 0009704
2018-06-06 07:56 peksa Note Added: 0009705
2018-06-06 08:50 henry Note Added: 0009706
2018-06-06 11:41 henry Note Added: 0009707
2018-06-06 12:21 peksa File Added: seulex_Imren.pdf
2018-06-06 12:21 peksa Note Added: 0009708
2018-06-06 12:36 henry Note Added: 0009709
2018-06-06 12:49 henry Note Added: 0009710
2018-06-07 07:49 peksa Note Added: 0009715
2018-06-07 08:04 henry Assigned To => henry
2018-06-07 08:04 henry Status new => resolved
2018-06-07 08:04 henry Resolution open => fixed
2018-06-07 08:04 henry Fixed in Version => 5.x
2018-06-07 08:04 henry Note Added: 0009716
2018-06-11 11:16 henry Status resolved => feedback
2018-06-11 11:16 henry Resolution fixed => reopened
2018-06-11 11:16 henry Note Added: 0009737
2018-06-12 14:45 peksa Note Added: 0009742
2018-06-12 14:45 peksa Status feedback => assigned
2018-06-22 12:08 henry Note Added: 0009811
2018-06-26 21:02 henry Status assigned => closed
2018-06-26 21:02 henry Resolution reopened => suspended
2018-06-26 21:02 henry Fixed in Version 5.x => dev
2018-06-26 21:02 henry Note Added: 0009833