Solver does the integration without calling the derivative callback function

Question

I have a python code (example from Cantera.org) that uses scipy.integrate.ode to solve a system of ODE. The code works fine and the results are reasnoable. However, I noticed something about the ode solver that does not make sense to me.

I have a put a print function inside (print("t inside ODE function", t)) the function the calculates the derivative vector (__call__(self, t, y)), and outside that function in the while loop (print("t outside ODE function", solver.t);). I expect that inside print has to be called when the solver does the time integration, and then the outside print is called. In other words, two "t outside ODE function" cannot appear right after another without "t inside ODE function" in between. However, this occurs in some of iterations in the while loop, which mean the solver does the integration without calculating the derivatives.

I am wondering how this is possible

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
        print("t inside ODE function", t)
        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('gri30.yaml')

# Initial condition
P = ct.one_atm
gas.TPX = 1001, P, 'H2:2,O2:1,N2:4'
y0 = np.hstack((gas.T, gas.Y))

# Set up objects representing the ODE and the solver
ode = ReactorOde(gas)
solver = scipy.integrate.ode(ode)
solver.set_integrator('vode', method='bdf', with_jacobian=True)
solver.set_initial_value(y0, 0.0)

# Integrate the equations, keeping T(t) and Y(k,t)
t_end = 1e-3
states = ct.SolutionArray(gas, 1, extra={'t': [0.0]})
dt = 1e-5
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("t outside ODE function", solver.t);
    print("\n")
    

Answer 1

The solver has an adaptive step size. Which means that it proceeds in internal steps that are adapted to the given error tolerances. In the segment from one step point to the next, the solution values get interpolated. Thus it can happen that a sequence of the external steps of the time loop falls into the same internal segment. If you set the error tolerances to smaller levels as the default ones, it can happen that the situation reverses, that several internal steps are required per external value request.