View Issue Details
ID | Project | Category | View Status | Date Submitted | Last Update |
---|---|---|---|---|---|
0002972 | OpenFOAM | Bug | public | 2018-06-05 11:07 | 2018-06-26 21:02 |
Reporter | peksa | Assigned To | henry | ||
Priority | low | Severity | minor | Reproducibility | always |
Status | closed | Resolution | suspended | ||
Platform | GNU/Linux | OS | Other | OS Version | (please specify) |
Fixed in Version | dev | ||||
Summary | 0002972: Stability / monotonicity check issue in the Seulex ODE solver | ||||
Description | In 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 Reproduce | I will provide an example of the issue later if needed. | ||||
Tags | ODE, seulex | ||||
|
It is not clear from your description what change you are proposing, could you clarify? |
|
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. |
|
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. |
|
> 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? |
|
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. |
|
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. |
|
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 |
|
I can confirm that your proposed implementation works for my cases as well. |
|
Resolved in OpenFOAM-5.x by commit 7f7d351b741bf6406366a043cac98de56d2d44dd Resolved in OpenFOAM-dev by commit 873d2c9c23a617432f96ca5cb9ffe6c02661e9c1 |
|
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. |
|
I will look into this tomorrow when I get back to my computer. |
|
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? |
|
Pending feedback from reporter |
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 |