function to return jacobian matrix

[Christopher Mauney]

Hello,

I know it's been asked before but I cannot find anything that suits what I'm looking for.

I want to use a another ODE integrator to integrate kinetics equations (similiar to the custom.py). However, my system will be very stiff and I would like to find the Jacobian matrix of the system.

First, I'm pretty sure this capability already exists in Cantera when it is internally integrating, so I am not sure if there is an interface I'm missing from the docs. (If so, or if not, this should be made clear).

Lacking this functionallity, I've considered:

1.) numerically calculating it: i.e, manually doing finite-differencing.

2.) some outside library (pyjac, numdifftools) 

As for (1.) its really tedious and inefficient. For (2.) the interfaces are poorly documented.

So, let's imagine I have a reaction network and some state:

gas = ct.Solution('h2o2.cti)

gas.TDX = 100.0, 100.0, 'H2:1, O2:1'

is there some call i could make to get the jacobian like:

f = gas.get_net_production_rates

J = gas.get_jacobian

Thanks

-Chris

[Nick Curtis]

Hi Chris,

I'm the lead developer of pyJac, and I'd love any input on how we can improve our documentation.

FWIW, we have a complete example of how to couple pyJac & cantera here.

As far as I know, it's not currently possible to extract the FD Jacobian from CVODE (the underlying ODE library that calculates the Jacobian for Cantera).

You could conceivably engineer this in, but it would require a lot of fiddling with the integrator interface (as well as further exposing it up the chain of the Reactor -> ReactorNet -> python interface).

I believe it's possible to extract the Jacobian for 1D flames, but I can't find it in the python documentation at the moment

Nick

[Ray Speth]

Chris,

You might be surprised, but where Cantera needs Jacobians internally, they are computed by finite difference. Making these internally-used Jacobians available would be a bit difficult.

In the case of reactor networks, the Jacobian is evaluated internally by CVODES, and I don’t think that CVODES provides an interface for accessing this matrix. In fact, I don’t think it really stores it for any length of time — I expect that it justs evaluate the Jacobian and then immediately factorizes it in-place.

In the case of the 1D solver, the Jacobian matrix is available via the C++ interface, via the jacobian() method of class OneDim (but again, the matrix is often in its factorized form). However, making this available in Python is non-trivial since I don’t think Numpy has a data structure equivalent to the BandMatrix class.

Regards,
Ray

On Monday, February 12, 2018 at 3:36:10 PM UTC-5, Christopher Mauney wrote:

[John D]

Has anyone done a one-on-one comparison of the Jacobian matrix computed by numdifftools and pyJac and compared the matrix elements?  If so, could the code segment to do so be shared with the user group?

Thanks

John

[solal....@gmail.com]

I have experience with both methods: finite difference and pyjac. Both yield approximately the same results but pyjac was shown to be much faster. I'll try to look for the code and share it here if/when I find it 

[John D]

Thanks.  That will be very helpful.  I tried using numdifftools with the gas.net_production_rates*gas.molecular_weights/gas.density to get dYi/dYj terms of the Jacobian and also followed the pyjac instructions to get  the Jacobian matrix (without the last species).  I see that all the Jacobian matrix elements with pyjac are 2x that produced by numdifftools (for all time-steps).  Any help from the community with more experience would be great. 

John

[Kyle Niemeyer]

Hi John,

I'm not familiar with numdifftools, but in our original study on pyJac, we actually found that it was *extremely* difficult to get high-fidelity Jacobians, at least for the chemical kinetics/reactor systems of interest, using finite differences (paper on arXiv, published version). We actually had to switch to an autodifferentiation tool, Adept, to verify the correctness of the Jacobian matrices.

Long story short: while there could certainly be bugs in pyJac, and we haven't been able to maintain/keep it updated, I would trust the results from it over finite differences, to be honest. Finite differencing is fine for approximate Jacobians, which is all you really need for implicit ODE integration (and you can get away with *very* approximate).

[John D]

Thanks, Kyle.  Here is the link to numdifftools

https://numdifftools.readthedocs.io/en/latest/

The test cases I tried with known analytical derivatives provide accurate results.

The temp profiles I get using pyJac with a time-integrator and the standard constant pressure reactor using cantera are the same.  So, I don't understand the 2x factor between the two Jacobians computed using pyJac and numdifftools.  I must be bungling something but can't figure it out.

Thanks

[Ingmar Schoegl]

Hi John,

One thing to keep in mind for Jacobians is that the equation of state puts constraints on how properties can change during a numerical perturbation. Depending on how you implement an analytical Jacobian, you can see results that are quite a bit different. Fwiw, Cantera's current development version (2.6.0a4) is introducing support for semi-analytical derivatives of various rates of progress and species creation rates. While these terms form building blocks for Jacobians, the actual implementation will depend on what type of reactor or reactor network you want to investigate (constant volume/pressure, or other).

-ingmar-

[John D]

Thanks, Ingmar.   I was looking at a regular ideal gas EOS with a pressure of 1 atm and 1200K initial temp.  Nothing fancy.  I started this exercise just because I was curious to see if I could take wdot, rho and molecular weights along Yi from Cantera and compute numerical Jacobians and see how well it compares with PyJac.  But after quite a bit of effort, I have not gotten very far.  I thought I would check with the Cantera community to see if someone more knowledgeable than me has done something like it.