sensBVP/examples/JacobianBasedSensitivity/jacVarSens.py

"""
Constant-pressure, adiabatic kinetics simulation with sensitivity analysis
"""

import sys
import time
import numpy as np

start = time.time()
import cantera as ct

gas = ct.Solution('./mech-FFCM1.cti')
temp = 900.0
pres = 10.0*ct.one_atm

gas.TPX = temp, pres, 'CH4:0.0499002,O2:0.199601,N2:0.750499'
r = ct.IdealGasConstPressureReactor(gas, name='R1')
sim = ct.ReactorNet([r])

# enable sensitivity with respect to all reactions
for i in range(0,gas.n_reactions):
    r.add_sensitivity_reaction(i)

# set the tolerances for the solution and for the sensitivity coefficients
sim.rtol = 1.0e-6
sim.atol = 1.0e-12
sim.rtol_sensitivity = 1.0e-6
sim.atol_sensitivity = 1.0e-12

#states = ct.SolutionArray(gas, extra=['t','s2'])
ignitionStateFound=False
Told=temp
TIgn=1.115367e+03 #K
tEnd=1.0  #s
tauIgn=0.0
nPts=1000
dt=tEnd/(float(nPts)-1.0)
sens=np.zeros(gas.n_reactions)
out=open("ignitionSensitivities.dat","w")

for t in np.arange(0, tEnd, dt):

    TOld=r.T
    sim.advance(t)
    TCurrent=r.T
    for i in range(0,gas.n_reactions):
        sens[i] = sim.sensitivity('temperature', i)

    if(ignitionStateFound==False):
        if(r.T>=TIgn):
            print("Ignition state found!\n")
            print("T=%15.6e K\n"%(TCurrent))
            print("tau=%15.6e s\n"%(t))
            ignitionStateFound=True
            tauIgn=t
            dTdtau=(TCurrent-TOld)/dt
            for i in range(0,gas.n_reactions):
                sens[i]=sens[i]*(-1.0e0/dTdtau)*(TCurrent/tauIgn)
                out.write("%15d\t%15.6e\n"%(i,sens[i]))
            break

out.close()
end = time.time()
print("Elapsed time: %15.6e s\n"%(end-start))