How to calculate more time in custom.py in cantera?

Question

I want to calculate auto-ignition delay time by using cantera in python linux. Here, I only change the t_end = 40 ， dt = 1，and gas.TPX = 801, P, 'H2:0.1667,O2:0.0833,N2:0.75' in custom.py example from Official Website. But, the result only calculate 20s. Anyone know how about this? my code:

import cantera as ct
import numpy as np
import scipy.integrate


class ReactorOde:
    def __init__(self, gas):
        # Parameters of the ODE system and auxiliary data are stored in the
        # ReactorOde object.
        self.gas = gas
        self.P = gas.P

    def __call__(self, t, y):
        """the ODE function, y' = f(t,y) """

        # State vector is [T, Y_1, Y_2, ... Y_K]
        self.gas.set_unnormalized_mass_fractions(y[1:])
        self.gas.TP = y[0], self.P
        rho = self.gas.density

        wdot = self.gas.net_production_rates
        dTdt = - (np.dot(self.gas.partial_molar_enthalpies, wdot) /
                  (rho * self.gas.cp))
        dYdt = wdot * self.gas.molecular_weights / rho

        return np.hstack((dTdt, dYdt))
gas = ct.Solution('h2o2.yaml')
P = ct.one_atm
gas.TPX = 801, P, 'H2:0.1667,O2:0.0833,N2:0.75'
y0 = np.hstack((gas.T, gas.Y))
ode = ReactorOde(gas)
solver = scipy.integrate.ode(ode)
solver.set_integrator('vode', method='bdf', with_jacobian=True)
solver.set_initial_value(y0, 0.0)
t_end = 40
states = ct.SolutionArray(gas, 1, extra={'t': [0.0]})
dt = 1
while solver.successful() and solver.t < t_end:
    solver.integrate(solver.t + dt)
    gas.TPY = solver.y[0], P, solver.y[1:]
    states.append(gas.state, t=solver.t)
print(states.t)

And the output is

[ 0.          1.          2.          3.          4.          5.
  6.          7.          8.          9.         10.         11.
 12.         13.         14.         15.         16.         17.
 18.         19.         19.45388009]

And the Fig about OH and T is here.

Answer 1

The loop terminates because the other condition, solver.successful() is not satisfied. The scipy.integrate.ode object also has a get_return_code() method that provides some hints about the problem. In your example, it returns -1, which, for the vode solver means:

Excess work done on this call. (Perhaps wrong MF.)

In this case, the "excess work" is because the intervals at which you're requesting output (1 second) are much longer than the time scales of the ODE that is being solved, which means that the solver has taken so many internal timesteps without reaching the next time that it assumes something has gone wrong.

Instead of specifying the output intervals explicitly, for stiff problems such as this one, you're better off just taking the solution at the times used by the integrator. This can be done by specifying the step=True argument to the integrate function, so it will just take one step toward the specified time. In this case, your integration loop would just look like:

while solver.successful() and solver.t < t_end:
    solver.integrate(solver.t + dt, step=True)
    gas.TPY = solver.y[0], P, solver.y[1:]
    states.append(gas.state, t=solver.t)

With this change, I see that the end time exceeds the specified t_end, as expected. You can also see that integrating this system requires time steps much smaller than the 1 second output interval you've specified by calculating the median step size:

>>> np.median(np.diff(states.t))
3.229414249972251e-05