Improve reactor network sensitivity analysis #55
Description
Abstract
Sensitivity analysis for reactor networks in Cantera does not seem to be as efficient or robust as it could be. I think this has something to do with how we are setting up the problem in CVODES, but what precisely the problem is is not obvious.
Motivation
Calculating sensitivities by solving the linear ODEs for the sensitivity coefficients should be a reliable and reasonably efficient method for computing sensitivities, compared to brute force, finite difference approaches. Currently, this is not the case.
Some examples demonstrating the current behavior, using the script appended at the end of this post (apparently you can't attach .py files). While this example calculates sensitivities with respect to enthalpy of formation, the same behavior can be observed for sensitivities with respect to reaction rate coefficients.
Agreement between CVODES and a finite difference approach are not great with default tolerances:
>>> idt = IgnitionDelayTest('h2o2.yaml')
>>> idt.test_ignition_delay_sensitivity()
Using CVODES sensitivities: 0.18 s
Using finite differences: 0.13 s
H2 -1.836805904e-08 -1.844576295e-08 -4.2215e-03
H 1.813382380e-08 1.817736038e-08 2.3980e-03
O2 -1.667790174e-09 -1.665175842e-09 -1.5688e-03
H2O2 -1.222888943e-09 -1.264490435e-09 -3.3450e-02
H2O 2.389793919e-10 2.387656571e-10 8.9477e-04
OH 3.603123634e-08 3.597439077e-08 1.5789e-03
HO2 -1.519046848e-08 -1.516084427e-08 -1.9521e-03The third column is the relative error in the ignition delay sensitivity to the enthalpy of formation of each species. Also note that even in this case, the finite difference approach is faster.
Tightening the sensitivity relative tolerance helps a little in this example (the default rtol_sens is 1e-4):
>>> idt = IgnitionDelayTest('h2o2.yaml', rtol_sens=1e-5)
>>> idt.test_ignition_delay_sensitivity()
Cantera version 2.5.0a4, git commit f41fdaf72
Using CVODES sensitivities: 0.20 s
Using finite differences: 0.13 s
H2 -1.845912950e-08 -1.844544233e-08 -7.4176e-04
H 1.815792709e-08 1.817721683e-08 1.0618e-03
O2 -1.665357045e-09 -1.665319392e-09 -2.2610e-05
H2O2 -1.261244315e-09 -1.264633985e-09 -2.6840e-03
H2O 2.389537735e-10 2.386221069e-10 1.3890e-03
OH 3.596316078e-08 3.597424722e-08 3.0822e-04
HO2 -1.515596387e-08 -1.516098782e-08 -3.3143e-04But tightening the tolerance further results in integrator failures. The effect of the sensitivity absolute tolerance is similar:
>>> idt = IgnitionDelayTest('h2o2.yaml', atol_sens=1e-9)
>>> idt.test_ignition_delay_sensitivity()
Using CVODES sensitivities: 0.18 s
Using finite differences: 0.13 s
H2 -1.845421440e-08 -1.844576295e-08 -4.5807e-04
H 1.817183559e-08 1.817736038e-08 3.0398e-04
O2 -1.665086245e-09 -1.665175842e-09 -5.3807e-05
H2O2 -1.262426139e-09 -1.264490435e-09 -1.6338e-03
H2O 2.389580004e-10 2.387656571e-10 8.0525e-04
OH 3.594692904e-08 3.597439077e-08 7.6366e-04
HO2 -1.516650911e-08 -1.516084427e-08 -3.7358e-04Further tightening either tolerance, or both together, tends to lead to never-ending integration, or integrator failures, e.g.:
>>> idt = IgnitionDelayTest('h2o2.yaml', rtol_sens=1e-5, atol_sens=1e-9)
>>> idt.test_ignition_delay_sensitivity()
Furthermore, sensitivity run time scales linearly with number of species. Sensitivities for 7 species:
>>> idt = IgnitionDelayTest('gri30.yaml', atol_sens=1e-5)
>>> idt.fuel = 'H2: 1.0, CH4: 1.0'
>>> idt.test_ignition_delay_sensitivity()
Using CVODES sensitivities: 0.92 s
Using finite differences: 0.79 sSensitivities for 53 species:
>>> idt.test_ignition_delay_sensitivity(species=idt.gas.species_names)
Using CVODES sensitivities: 7.92 s
Using finite differences: 5.42 sPossible Solutions
A few things come to mind that may be worth investigating:
-
With the current setup, it seems that CVODES is forced to evaluate derivatives of the residual function with respect to the perturbed reaction rates one at a time. Is there a way to avoid this? These should be fairly simple analytical derivatives, at least for the reaction rate sensitivities. Perhaps less so in the case of enthalpies, but even then, only a handful of reaction rates / state variables should depend on any one enthalpy of formation, so we could still reduce the cost of this evaluation.
-
Would this improve if we were using sparse/iterative methods, even without analytical derivatives? Could we expect this to be resolved efficiently as an extension to Addition of Jacobian Preconditioning Functionality to Cantera. #47?
-
Would having derivatives calculated using automatic differentiation help? This is a fairly big project, but one that would have a lot of benefits in other areas as well.
Example Script
This example is adapted from the Python test test_reactor.TestReactorSensitivities.test_ignition_delay_sensitivity.
import cantera as ct
import numpy as np
from timeit import default_timer as timer
print(f'Cantera version {ct.__version__}, git commit {ct.__git_commit__}')
class IgnitionDelayTest:
def __init__(self, mech, rtol=1e-9, atol=1e-15, rtol_sens=1e-4, atol_sens=1e-6):
self.gas = ct.Solution(mech)
self.t_end = 0.6
# Default values correspond to Cantera defaults
self.rtol = rtol
self.atol = atol
self.rtol_sensitivity = rtol_sens
self.atol_sensitivity = atol_sens
self.phi = 0.4
self.fuel = 'H2: 1.0'
self.oxidizer = 'O2:1.0, AR:4.0'
def setup_ignition_delay(self):
self.gas.TP = 900, 5*ct.one_atm
self.gas.set_equivalence_ratio(self.phi, self.fuel, self.oxidizer)
self.reactor = ct.IdealGasReactor(self.gas)
self.net = ct.ReactorNet([self.reactor])
self.net.rtol = self.rtol
self.net.atol = self.atol
self.net.rtol_sensitivity = self.rtol_sensitivity
self.net.atol_sensitivity = self.atol_sensitivity
def calc_tig(self, species, dH):
self.setup_ignition_delay()
S = self.gas.species(species)
st_orig = S.thermo
coeffs = st_orig.coeffs
coeffs[[6, 13]] += dH / ct.gas_constant
s_new = ct.NasaPoly2(st_orig.min_temp, st_orig.max_temp, st_orig.reference_pressure, coeffs)
S.thermo = s_new
self.gas.modify_species(self.gas.species_index(species), S)
t = []
T = []
while self.net.time < self.t_end:
t.append(self.net.time)
T.append(self.reactor.thermo.T)
self.net.step()
T = np.array(T)
t = np.array(t)
S.thermo = st_orig
self.gas.modify_species(self.gas.species_index(species), S)
return (t[-1]*T[-1] - np.trapz(T,t)) / (T[-1] - T[0])
def calc_dtdh(self, species):
self.setup_ignition_delay()
for s in species:
self.reactor.add_sensitivity_species_enthalpy(s)
t = [0.0]
T = [self.reactor.T]
S = [[0.0]*len(species)]
iTemp = self.reactor.component_index('temperature')
while self.net.time < self.t_end:
self.net.step()
t.append(self.net.time)
T.append(self.reactor.thermo.T)
S.append(self.net.sensitivities()[iTemp])
self.T = T = np.array(T)
self.t = t = np.array(t)
self.S = S = np.array(S)
To = T[0]
Tf = T[-1]
tig = (t[-1]*Tf - np.trapz(T,t))/(Tf-To)
dtdp = ((t[-1] - tig)*S[-1,:]*Tf - np.trapz(S*T[:,None], t, axis=0))/(Tf-To)
return dtdp
def test_ignition_delay_sensitivity(self, species=('H2', 'H', 'O2', 'H2O2', 'H2O', 'OH', 'HO2')):
clock1 = timer()
dtigdh_cvodes = self.calc_dtdh(species)
clock2 = timer()
print(f'Using CVODES sensitivities: {clock2 - clock1:.2f} s')
clock1 = timer()
tig0 = self.calc_tig('H2', 0)
dtigdh = np.zeros(len(species))
dH = 1e4
for i,s in enumerate(species):
dtigdh[i] = (self.calc_tig(s, dH) - tig0) / dH
clock2 = timer()
print(f'Using finite differences: {clock2 - clock1:.2f} s')
for i,s in enumerate(species):
print('{:4s} {: .9e} {: .9e} {: .4e}'.format(
s,
dtigdh_cvodes[i],
dtigdh[i],
2*abs(dtigdh_cvodes[i] - dtigdh[i]) / (dtigdh_cvodes[i] + dtigdh[i])))
if __name__ == '__main__':
idt = IgnitionDelayTest('gri30.yaml', atol_sens=1e-5)
idt.fuel = 'H2: 1.0, CH4: 1.0'
idt.test_ignition_delay_sensitivity(species=idt.gas.species_names)