Problem description

Depending on the output times given to Reactor.advance, specified limits on the change in one or more state variables can be exceeded by a wide margin. In some cases, this seems to stem from the relevant variable having a near-zero derivative instantaneously, permitting a large interval to be passed on to the internal advance method, which then takes many time steps over which the variable changes by more than the expected amount. However, there are other cases where it seems like the instantaneous rate of change should predict the need for a shorter time step.

I originally encountered this due to surprising failures for TestReactor.test_advance_with_limits in #1452, where it seems that the fact that test passes at all is kind of a fluke.

Given the intention for this feature to be a robust way of getting nice output even for inexperienced users (see discussion in #628/#629), I think some improvement here is warranted.

Steps to reproduce

The following is based on TestReactor.test_advance_with_limits, with a few modifications.

import cantera as ct
import matplotlib.pyplot as plt
gas = ct.Solution('gri30.yaml')

def integrate(limit=None):
    P0 = 10 * ct.one_atm
    T0 = 990
    X0 = 'CH4:1.0, O2:1.0, AR:5.0'
    gas.TPX = T0, P0, X0
    r1 = ct.IdealGasReactor(gas)
    net = ct.ReactorNet([r1])
    ix = net.global_component_index('CH4', 0)
    r1.set_advance_limit('CH4', limit)

    tEnd = 1.0
    tStep = 0.1
    nSteps = 0

    states = ct.SolutionArray(gas, extra=['t', 'steps'])
    t = tStep
    yp = net.get_state()[ix]
    steps_p = net.solver_stats['steps']
    dy_max = 0
    while t < tEnd:
        steps = net.solver_stats['steps']
        states.append(TPX=r1.thermo.TPX, t=net.time, steps=(steps-steps_p))
        steps_p = steps
        t_curr = net.advance(t)
        y = net.get_state()[ix]
        dy_max = max(dy_max, abs(y-yp))
        yp = y
        nSteps += 1
        if t_curr >= t:
            t += tStep

    print(f'case: {limit=}, {dy_max=}, {nSteps=}')
    return nSteps, states

n_baseline, states_baseline = integrate()
n_advance_coarse, states_coarse = integrate(0.005)
n_advance_fine, states_fine = integrate(0.001)

f, ax = plt.subplots(figsize=(5,3))
ax.plot(states_fine.t, states_fine('CH4').Y, '.-', label='fine')
ax.plot(states_coarse.t, states_coarse('CH4').Y, '.-', label='coarse')
ax.legend();

Behavior

case: limit=None, dy_max=0.06457704195575183, nSteps=10
case: limit=0.005, dy_max=0.0645770419565304, nSteps=10
case: limit=0.001, dy_max=0.014385584336446846, nSteps=90

As you can see, in the second case, the desired limit is exceeded by more than a factor of 10, and the output jumps over the ignition event entirely. Looking at states_coarse.steps, you can see that integrator takes over 2000 time steps between the second and third output points.

With a finer limit, the behavior is better, but this isn't even a guaranteed path to success -- if tStep is increased to 0.2 s, then both values of the advance limit miss everything. Also, even here, there is a strange gap toward the end of the ignition event where the fuel mass fraction drops by quite a bit with no intervening output.

System information

• Cantera version: 2.6.0 or current main branch at 399e1cb
• Python 3.10 used for the example above, but affects all versions