Potential issue with cython wrapper

[michael-a-hansen]

Hello,

I've just used pyJac's pyWrap utility to write C code that I can call from python. I'm getting nans when I evaluate py_dydt and py_eval_jacobian, but not from py_eval_conc, py_eval_rxn_rates, py_get_rxn_pres_mod, or py_eval_spec_rates. See details below.

Build procedure

pyJac call: python3 -m pyjac -i h2-burke.xml -l c pyWrap call: python3 -m pyjac.pywrap --source_dir out -l c

While pyWrap completes without errors, I get a deprecated numpy warning during compilation and the following warning during linkage. I'm not sure but it seems like it might be important.

ld: warning: ignoring file /Users/mike/Mechanisms/coh2-davis/build/temp.macosx-10.6-intel-3.6/libc_pyjac.a,
file was built for archive which is not the architecture being linked (i386):
/Users/mike/Mechanisms/coh2-davis/build/temp.macosx-10.6-intel-3.6/libc_pyjac.a

Here's the input file (XML generated by cantera): h2-burke.xml.zip

Sample code (requires cantera, numpy, and generated pyjacob)

The snippet below evaluates some pyJac functions and prints the results, demonstrating (on my machine) that everything looks okay except for dYdt and Jflat, which are full of nans.

import pyjacob as pyjacob
import cantera as ct
import numpy as np

T = 300
p = 101325
gas = ct.Solution('h2-burke.xml')
nspec = gas.n_species
Y = 'H2:2, O2:1, N2:3.73'
Y = np.ones(nspec)
gas.TPY = T, p, Y

Y = gas.Y
mw = gas.mean_molecular_weight
rho = gas.density_mass

conc = np.zeros_like(Y)
dYdt = np.zeros(nspec + 1)
molar_sp_rates = np.zeros(nspec)
Jflat = np.zeros((nspec + 1) * (nspec + 1))

fwd_rates = np.zeros(gas.n_reactions)
rev_rates = np.zeros(gas.n_reactions)
pres_mod = np.zeros(gas.n_reactions)

pyjacob.py_eval_conc(T, p, Y, mw, rho, conc)
pyjacob.py_eval_rxn_rates(T, p, conc, fwd_rates, rev_rates)
pyjacob.py_eval_rxn_rates(T, p, conc, fwd_rates, rev_rates)
pyjacob.py_get_rxn_pres_mod(T, p, conc, pres_mod)
pyjacob.py_eval_spec_rates(fwd_rates, rev_rates, pres_mod, molar_sp_rates)

print(T, p, mw, rho)
print(Y)
print(conc)
print(fwd_rates)
print(rev_rates)
print(pres_mod)
print(molar_sp_rates)

pyjacob.py_dydt(T, p, Y, dYdt)
pyjacob.py_eval_jacobian(T, p, Y, Jflat)
print(dYdt)
print(Jflat)

Software state

• Mac OS X 10.12.3, XCode 8.2.1
• pyJac updated today
• python 3.6 with numpy, cython, cantera all updated

[skyreflectedinmirrors]

@michael-a-hansen ok, the issue here is that the t argument in py_dydt and py_eval_jacobian actually represents time, not temperature. The state vector (for the C module) should consist of:

[T, Y0, Y1.... YNsp]

Hence you should get the correct results by calling this as:

Y = gas.Y
Y = np.concatenate(([T], Y))
dYdt = np.zeros_like(Y)
Jflat = np.zeros((nspec) * (nspec))
pyjacob.py_dydt(0, p, Y, dYdt)
pyjacob.py_eval_jacobian(0, p, Y, Jflat)

Also note that due to the formulation of our system, the Jacobian is actually Nsp x Nsp, i.e. 1 entry for Temperature and Nspecies - 1 entries for the species mass fractions (using the strict mass conservation formulation)

The reason for all this confusion is that typically in ODE solvers, the RHS / Jacobian function accept time as an argument in case the equations are dependent on it directly. To conform to CVODE, etc. it is also included here.

@kyleniemeyer didn't we have someone working on documentation for the python wrapper? It looks like an example would be particularly useful 👍

edit: also, we might want to just drop the 'time' argument, although that may hamper coupling pyJac to a python based integrator package

[michael-a-hansen]

Terrific! Thank you. It is working as expected now :)

A note: it's very confusing that Y is TY in the higher-level functions (dydt and such) and Y in the lower-level methods (eval_conc and such). Additionally, when looking at the auto-generated code, I saw a line like the following, which I was ready to submit as a massive bug report before I realized Y is TY and so y[1] is actually the first mass fraction (first element of Y, second of TY)! I think simply replacing Y with TY where appropriate would help new developers and users immensely.

cp[0] * y[1] + cp[1] * y[2] + cp[2] * y[3] + cp[3] * y[4] ...