Round of error in zeroD reactor simulation wih Python interface

[ali shahanaghi]

Dear all,

I recently observed the below deviation in the results of Cantera (python interface) zeroD reactor simulations.

So I ran the same code (reactor1.py) using different versions of Cantera and on different machines. I observed different floating point values in order of 1e-6 for temperature and 1e-4 for internal energy.

I was wondering if this is a known issue and if there is any way to avoid that. 

My simulations concern using ZND model (SDToolbox) in which the chemical length/time scales are in order of 1e-6. Therefore such differences sometimes make convergence issues especially when I use a detailed mechanism.

Best Regards,

Ali Shahanaghi

[Bryan Weber]

Hi Ali,

The controls for the tolerances are atol and rtol on the Reactor Network or set_steady_tolerances and set_transient_tolerances for 1-D problems.. Roughly, rtol (relative tolerance) gives you the number of significant figures in the result for large numbers, whereas atol (absolute tolerance) takes that role for small numbers. This page has a pretty good description of the difference for a generic Matlab ODE solver. The same concepts apply to the CVODE solver used by Cantera.

In your case, you may want to decrease both those numbers, keeping in mind the limitations of floating point numbers, see: https://wiki.sei.cmu.edu/confluence/display/c/FLP00-C.+Understand+the+limitations+of+floating-point+numbers and https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html As you can also see from those links (and your results) the exact answer also depends on how Cantera is compiled (that is, the specific compiler and compiler optimizations used).

Finally, I doubt (without evidence except experience) that differences on the order that you show in your screenshot will affect integrated/global responses substantially. Aside from the fact that global responses are generally not extremely sensitive to the tolerances, you're very likely to have significant errors in your detailed chemistry that actually control the accuracy of your solution relative to experimental results. That said, if you're concerned about chemical time scales on a particular order in a transient simulation, you should adjust atol and rtol so that CVODE chooses an appropriate time step... if you're concerned about length scales in a 1-D simulation, you should choose appropriate grid refinement criteria as well as the tolerances, so that the solver resolves the length scale you're interested in.

Hope that helps,

Bryan

[ali shahanaghi]

Thank you very much Bryan.

I initially found the default rel/abs tolerances set in CVodesIntegrator (1e-9/1e-15) to be low enough, but, I will test lower values as well.

Below there is a screenshot from Cantera 2.5.1 using Anaconda precompiled installation. I set gas.TP and then print gas.TP again.

As you can see the pressure value has changed after initialization.

Could it be from the compiler floating point limitation as you mentioned? or is there any operation in the background on the pressure that changes its value?

Thanks again,

Ali

[Bryan Weber]

Hi Ali,

The natural state variables for an ideal gas in Cantera are the temperature and the density. So when you set the state with the temperature and pressure, Cantera has to do some arithmetic to store the appropriate value of the density. You can see this here: https://cantera.org/documentation/docs-2.6/doxygen/html/de/dbd/IdealGasPhase_8h_source.html#l00371 where the pressure method returns the calculated value of the pressure based on the ideal gas equation, and setting the pressure is implemented by setting the density based on the ideal gas equation.

So the difference that you see is likely the accumulation of a bit of floating point error in those calculations, plus some floating point error in representing the initial value of the pressure.

Just to say again, if your results are so sensitive to these kinds of differences, then it seems to me you're unlikely to be able to use any numerical schemes to solve your equations. Also, the error introduced into the problem by incorrect/inaccurate rate constants is almost certainly going to be far larger than the numerical error from floating point rounding, relative to the physical "truth" of the system (if you could conduct an physical experiment).

Hope that helps,

Bryan

[ali shahanaghi]

Thank you, Bryan. 

Very helpful instructions. Indeed, the deviation between provided rate constants and experiments is a known issue in this type of simulation.

My confusion was mainly about running the code in different ways (e.g. sometimes in a loop for reaction sensitivity analysis or on different machines) 

and getting results that were different (or diverging).

Now with provided hints, I can find a workaround to reduce the sensitivity of my results.

Regards,

Ali