Stiff integration of the chemical source term

[Solal Amouyal]

I have a 3D reacting CFD code (Fortran) which uses Cantera for computing the thermodynamic variables (viscosity, heat capacity, etc...) and the kinetics (net production rates, ...).

However for the integration of the chemical source term in the species and energy equation, the integration itself is done using an external solver. I was wondering, would it possible (and efficient) to do the integration in Cantera itself? Looking around, I see that the library provides the Reactor class (not provided for the Fortran interface but that should not be a problem). Is that efficient enough to use as a stiff integrator in a CFD code where CPU time is critical or would you recommend sticking to standard solvers such as VODE or CVODE?

[Kyle Niemeyer]

Hello Solal,

The ReactorNet (https://www.cantera.org/docs/sphinx/html/cython/zerodim.html#reactor-networks) class actually uses the CVODE integrator from SUNDIALS to advance its state, so for most problems you can just rely on that.

(Also, in general I'd recommend using that over the older VODE, as it should offer better performance.)

Kyle

[Nick Curtis]

Solal,

This is a somewhat in-depth topic, but some general rules to guide you:

1)  VODE / CVODE has a large startup time.  When it is "restarted" (i.e., at the beginning of each ODE integration per-cell) it uses a small starting time-step and slowly builds up solution accuracy and convergence order.  This is a good choice to handle large amounts of stiffness, but may be overkill for small time-steps (e.g., for LES).  Other ways to speed-up CVODE include using a good pre-conditioner (see McNenly & Whitesides "Adaptive Preconditioning Strategies for Integrating Large Kinetic Mechanisms") but this is often problem/mechanism dependent.

As Kyle mentioned, CVODE is used internally by Cantera already

2) Extrapolation methods, e.g., SEULEX / SIBS are another viable option.  These have much lower start-up costs, but accuracy and efficiency may suffer for the tolerances you're going to be looking at for CFD (https://scicomp.stackexchange.com/questions/27178/bdf-vs-implicit-runge-kutta-time-stepping)

2)  Next up are so-called linearly implicit methods (i.e., they do not require newton iteration).  These have essentially zero start-up cost, but require an accurate Jacobian to maintain accuracy / order.  This means you should either use auto-differentiation or an analytical Jacobian e.g., pyJac (https://github.com/SLACKHA/pyJac)

Best,

Nick

[Solal Amouyal]

Kyle, Nick,

Thank you for the replies.

I was aware of the different methods mentioned, but not familiar with them, which is why I favored the use of a "black-box" solver which could take care of the specifics for me (choice of methods or preconditioning/Jacobian etc...).

However I was not aware of your package pyJac which I think will prove very useful! I am working on a high-order LES solver so the third option looks the most adequate for my use. I will look in the literature for a specific scheme but if you know of one which you works well for hydrocarbon combustion, I will be happy to hear it. 

Sincerely,

Solal 

[Nick Curtis]

Hi Solal,

The best scheme to use depends on your case, i.e., the stiffness of the mechanism / timestep of your solver, and unfortunately a lot of them don't have "standard" implementations outside of the Harrier and Wanner fortran codes.

If you are using OpenFOAM by chance, there are a bunch of Rosenbrock type solvers already implemented (orders 2-4, as well as some similar Rodas type solvers and SIBS / Seulex). Boost's ODEInt library has a 4th order rosenbrock solver implementation as well.

Best,

Nick

[Nick Curtis]

Also, we've developed an exponential rosenbrock-type integrators (w/ krylov subspace approximation) that might be useful for this kind of time-step.

See this paper for algorithm details, and here for the implementation study (note: we were particularly looking at GPU-based integrators w/ somewhat larger time-steps O(1e-6)).

These solvers can use pyJac by default.

We are planning on incorporating rosenbrock type solvers for our new release (i.e., once I get it linked up to our new version of pyjac)

Nick

[Solal Amouyal]

Hi Nick,

I am familiar with your (and Kyle's) work but I admit that I had the false impression the open-source software you were developing were solely dedicated to GPU computing. 

The installation of the accelerInt solver failed so I focused on Boost's solvers.

I wrote an example C++ script which integrates a H2-O2 mixture with Boost's ODEint library and the scripts generated by pyjac. I was able to get good results with an implicit RK method but the Rosenbrock integrator gets stuck or takes too long even with low tolerances. Compared with a Euler integration (my own), the RK scheme is still substantially slower; the mechanism used (9 species, 28 reactions) is probably too small to see the benefits of implicit methods. 

My main objective is to run my CFD code with DME combustion, which entails more complex mechanisms (39 species, 175 reactions). I will try to see if I get more acceleration which that mechanism.

I have attached my code if someone is also interested with this kind of work. Just change the cantera installation directory in the source/Makefile.

Sincerely, 

Solal

[Nick Curtis]

Hi Solal,

It's worth noting that one of the other users in our group reported similar problems to yours a few weeks ago -- slow integration / getting stuck w/ the linearly implicit methods, but not the implicit RK -- the solution in that case was to simply transpose the Jacobian, at which point the slowdown for (in that case) the exponential integrators.

There was an error in the pyJac documentation regarding the ordering of elements in the Jacobian (though, the website has since been updated), I have a a few pull requests regarding this I need to finish up.

The accelerInt documentation is also woefully incomplete (again, I've been meaning to fix this up but time is always an issue!), but if you post your error message (or open an issue on GH) I'd be happy to take a look.

I'll try to take a look at your case in the next few days, but I suspect it's a similar issue.

Best, 

Nick

[Nick Curtis]

Solal,

I decided to take a quick look at your code, the problem I noticed is that your Jacobian is sized (Nsp + 1)x(Nsp + 1), however the Jacobian in pyJac specifically excludes the last species (as all those derivatives are identically zero), and should be sized (Nsp)x(Nsp).

I've attached an updated version of the main.cpp file that actually runs, and gives about what you'd expect:

Rosenbrock4 results:
       - Number of steps:     3802
       - CPU time = 5.310e+00
              0                  2.9328e+03
              1                  5.5290e-03
              2                  1.2829e-03
              3                  7.1677e-03
              4                  2.2285e-02
              5                  2.6598e-02
              6                  1.9313e-01
              7                  1.2844e-05
              8                  7.9674e-07
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
RK-dopri5 results
       - Number of steps:   516661
       - CPU time = 2.974e+01
              0                  2.9328e+03
              1                  5.5290e-03
              2                  1.2829e-03
              3                  7.1677e-03
              4                  2.2285e-02
              5                  2.6598e-02
              6                  1.9313e-01
              7                  1.3240e-05
              8                  3.6862e-07
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Euler integration:
       - CPU time = 3.610e+00
              0                  2.9328e+03
              1                  5.5290e-03
              2                  1.2829e-03
              3                  7.1678e-03
              4                  2.2285e-02
              5                  2.6598e-02
              6                  1.9313e-01
              7                  1.2175e-05
              8                  1.3903e-06
              9                  7.4399e-01
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Cantera equilibrate:
              0                  2.9328e+03
              1                  5.5290e-03
              2                  1.2829e-03
              3                  7.1677e-03
              4                  2.2285e-02
              5                  2.6598e-02
              6                  1.9313e-01
              7                  1.2844e-05
              8                  7.9674e-07
              9                  7.4399e-01

I'd imagine for the stiffer DME case, the Euler integration would start to fall apart (or you would have to lower the time-step substantially)

Best,

Nick

[Kyle Niemeyer]

Solal-

I'll just briefly follow up by saying that, since this issue relates to ongoing efforts outside of Cantera by Nick and I, we'd be happy to work with you on this more one-on-one, if you want to move the conversation outside of the user group.

Kyle