[Python] Fix ATOL in advance_to_steady_state

[rwest]

The docstring says that if atol is undefined it should match the solver ATOL,
which makes sense to me. What the code in fact did was atol = self.rtol,
which I presume is a typo, so have changed to atol = self.atol.

Since RTOL is often ten or more orders of magnitude larger than ATOL,
this could make some significant differences to people trying to use this
method.

However, the changes could be bad: if now we can't ever reach the steady state,
or it takes much longer to reach steady state (already a complaint I hear),
folks may complain that this commit just made it a bunch slower and less stable.

(But at least they're now getting what they ask for, instead of just getting
the wrong answer quicker. ).

Changes proposed in this pull request

• Set atol = self.atol in advance_to_steady_state as described in the method's docstrings.

Checklist

• The pull request includes a clear description of this code change
• Commit messages have short titles and reference relevant issues
• Build passes (scons build & scons test) and unit tests address code coverage
• Style & formatting of contributed code follows contributing guidelines
• The pull request is ready for review

What I haven't yet done is show what effect this has, but opening the PR now for comments.

[ischoegl]

This indeed does look like a typo - I believe the suggested fix is appropriate.

[12Chao]

However, the changes could be bad: if now we can't ever reach the steady state,
or it takes much longer to reach steady state (already a complaint I hear),
folks may complain that this commit just made it a bunch slower and less stable.

This does indeed make some situations crash that didn't use to crash, I was running a methane partial oxidation on Pt model simulation with a chain of CSTR, the tolerances I used was rtol=1e-10 and atol=1e-20, and I got a CVODE error message after I corrected the atol.

CanteraError: 
*******************************************************************************
CanteraError thrown by CVodesIntegrator::step:
CVodes error encountered. Error code: -27

tout too close to t0 to start integration.
Components with largest weighted error estimates:
0: -1.7823798662745834e-05
2: 8.137304116665512e-06
6: 5.089364868104281e-06
7: -3.104365468736558e-06
23: 1.820722907327393e-11
14: 1.0850536718571933e-12
10: 2.9499388397008586e-13
16: 1.678002681128422e-13
13: 1.7548594757486803e-15
1: 1.3434606588734786e-15
*******************************************************************************

I guess it's expected because the actual atol makes it stiffer, but that's probably a good reason to justify that an improvement to the steady-state solver is needed. There is an issue #654 about this. Since I came across this too, I will take a look at #1021, see if I could help on that one.

[ischoegl]

@12Chao ... did you try passing a reduced atol as a parameter to advance_to_steady_state?

[12Chao]

No, I changed atol=self.rtol to atol=self.atol in Cantera and ran advance_to_steady_state().

[ischoegl]

I believe that using advance_to_steady_state(atol=??) should fix your issues, see documentation

[bryanwweber]

1e-20 is really small for atol, that represents roughly the number of significant figures... A standard double precision number has about 15 sig figs. So you don't need to make atol that small.

[12Chao]

1e-20 is really small for atol, that represents roughly the number of significant figures... A standard double precision number has about 15 sig figs. So you don't need to make atol that small.

Thanks, that's good to know, I guess I'll try to avoid using very small atol in the future simulations, but a 1e-20 atol worked fine with advance method, I replaced advance_to_steady_state to advance(sim.time + 10 * residual_time), and the simulation was finished smoothly. It makes me wonder why advance can integrate with smaller atol, is that because advance does not use step? could you please give me some hint on that?

[bryanwweber]

I would bet that the problem is that advance_to_steady_state is advancing to a much larger time than 10*residual_time (do you mean residence time?). The advance_to_steady_state uses some scaling method to try and determine when the simulation has reached a steady state, and keeps running advance until the criteria are satisfied. I suspect that the criteria aren't being satisfied, and CVODES can't keep integrating for such a long time with such tight tolerances. The criteria were developed for homogeneous problems, as far as I know, so they may not be tuned for heterogeneous problems. So the problem is not the atol per-se, but rather that the advance_to_steady_state algorithm isn't robust to this failure mode. You can check this by printing the simulation time when the failure occurs (probably using try/except is the easiest way), hopefully I'm not way off base here 😄