Jacobians of Kinetic Source Terms

[ischoegl]

I wanted to start a discussion on Jacobians, as there is ongoing activity (Cantera/cantera#1010), as well as some work-in-progress (#19) and external pertinent code (e.g. pyJac). As an aside, I have coded up analytical Jacobians myself in the past, although that code is not released. Further, there are obvious extensions that go beyond what I want to discuss here (adjoint sensitivity analysis, automatic differentiation), so I mainly want to focus on how to make the exact evaluation of Jacobians most useful to the Cantera community. Per recent comment by @kyleniemeyer, finite difference implementations are often insufficiently accurate, so I believe this would be an important feature.

To me, the main question is whether the evaluation of a Jacobian should be available as a method of Solution (i.e. tied into Kinetics.h, with different options, akin to equilibrate), or whether it would make sense to provide derivatives that allow users to assemble Jacobians (fyi, there's a plug for a temperature derivatives of reaction rates from my recent work here). One question here is how to ensure that the sparsity is available at various front ends (e.g. via scipy.sparse?). This obviously extends to Jacobians for entire reactor networks, but is likewise relevant to source terms within one-dimensional flames.

I'd appreciate if this discussion could take place on a public venue, and I'm tagging people who may want to chime in: @kyleniemeyer, @anthony-walker, @speth, @DavidAkinpelu and @mbkumar (there may be others).

PS: Adding some references:

• Analytical Jacobian for Cantera #19 - analytical Jacobians
• Addition of Jacobian Preconditioning Functionality to Cantera. #47 - preconditioning for Sundials solvers (Adding adaptive preconditioner functionality to Cantera cantera#1010 / formerly Adding adaptive preconditioner functionality to Cantera (draft) cantera#951)
• Use automatic differentiation to calculate derivatives #80 - automatic differentiation
• Provide Jacobian to Sundials solvers #86 - Jacobians for Sundials solvers
• Refactor Reaction Rate Evaluators #87 - Refactoring of reaction rate evaluators (also Refactor reaction rate evaluators cantera#995)
• Reconsider composition state variables for reactors, flames, and other systems #91 - Discussion of alternative composition state variables

[ischoegl]

Application Example - Jacobian of Solution 'rates'

As a starting point, a Kinetics-centered approach could work as follows (I am using Python to illustrate my point)

import cantera as ct
gas = ct.Solution(`some_mechanism.yaml`)

... # perform some manipulations on gas

jac = gas.jacobian(mole=True, heat_release=True, sparse=True) # placeholders; exact interface tbd

This would allow for access that roughly maps to existing properties gas.net_production_rates and gas.heat_production_rates, which are the standard kinetic source terms in governing equations, and have obvious partial derivatives with respect to species concentrations, temperature, and pressure. There are obviously some open questions (e.g. how to handle pressure?, stacking order?, separate Jacobians?), which remain to be discussed.

Anchoring this capacity within Kinetics would imho not restrict things for other expansions (e.g. Jacobians for reactor networks), as it would be straight-forward to loop through kinetics managers to populate the composite Jacobian matrix. Essentially, this would use the same approach as what is already done to populate the composite solution vector; reactor networks objects don't implement specifics of individual source terms, but just assemble the governing equation from parts. Caveat: there may be some chain rules based on the choice of conserved properties.

The point of this discussion is to ensure that any implementation would provide the most flexible interface for the Cantera community (as well as efficient calculations internally). I'd suggest to keep replies to comments application-specific, i.e. reply below with responses to this application, but start a new comment for the discussion of another application.

[ischoegl]

@kyleniemeyer and @anthony-walker

As we need Kinetics-based (sparse) Jacobians in our own work (for what we need, pyJac doesn't apply, and spyJac installation appears to be broken), I started to implement sparse analytic species derivatives in PR Cantera/cantera#1081 from scratch. At this point, the following are implemented for GasKinetics, and are accessible from Python via:

>>> gas.rop_species_derivatives(forward=True, reverse=True)
>>> gas.production_rate_species_derivatives(creation=True, destruction=True)

I still need to implement Jacobian terms accounting for third-body efficiencies (which is straight-forward) and concentration-dependent rate terms (which I haven't looked into yet). What is implemented is covered by unit tests, which check for equivalency with a numerical derivative.

[anthony-walker]

Thanks for letting me know about this @ischoegl. The profiling that I did on Cantera/cantera#1010 revealed that there was some performance impacts in the way that I am calculating the species derivatives with the maps. I have started working on re-implementation of the method for now but it seems you have done quite a bit here in that regard and I don't necessarily want to repeat your work. The work I am doing actually neglects third body efficiencies as one of it's simplifications so it may already be in a useful state for what I need.

[ischoegl]

Sure - we're interested in exact Jacobians, so this isn't quite where we need it. I'll likely add third bodies to that PR, but may leave the rest for later so this can be reviewed and merged quickly.

[ischoegl]

Cantera/cantera#1089 supersedes Cantera/cantera#1081. It is nearly feature complete at this point, but I expect some discussion as our own use case is not the only relevant context where Jacobians are relevant.

[ischoegl]

🎉 ... Cantera/cantera#1089 is now merged.