Reactor.advance limits can be exceeded greatly

[speth]

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

[ischoegl]

@speth ... the problem for the surprising behavior is that tStep is excessively large - as implemented, the adaptive refinement relies exclusively on predictions based on currently available derivatives. I agree that this can fail in some instances, where it may be necessary to recalculate intermediate results (which the current implementation seeks to avoid).

[speth]

tStep in the example above is actually smaller (relative to the ignition delay time) than it is in the current unit test. I think the only reason that the unit test sort of works is because the advance time step is precisely aligned with the IDT. If I take the conditions from the unit test, just reducing the end time so that things actually show up on a plot, these are the results I get:

case: limit=None, dy_max=0.005447321778656981, nSteps=9
case: limit=0.01, dy_max=0.005238599220587466, nSteps=13
case: limit=0.001, dy_max=0.0018037034355080684, nSteps=23

Full code for this example

gas = ct.Solution('h2o2.yaml')

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

    tEnd = 0.01
    tStep = 1e-3
    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.01)
n_advance_fine, states_fine = integrate(0.001)

While it does take a few more steps in the "coarse" case with the limit imposed, they're not really in the right place. And if you pick time steps that aren't aligned with the IDT, even if they are much smaller than the IDT, then the "fine" case also tends to struggle. For example, setting:

    tEnd = 0.01
    tStep = 1.1e-4
    nSteps = 0

Gives the results:

case: limit=None, dy_max=0.005515943151123514, nSteps=90
case: limit=0.01, dy_max=0.005365298679915316, nSteps=92
case: limit=0.001, dy_max=0.005365298679915316, nSteps=93

Here, you can see that the output is skipping over almost the entire ignition process despite a fair number of steps occurring before the event.

I don't think it's necessary to implement backtracking to fix this. In all of these cases, CVODES is taking several hundred timesteps internally, within one call to CVode with the CV_NORMAL flag. If the advance constraint was checked somewhat more frequently (probably not even every timestep), I think this could be resolved.

[ischoegl]

Looking at the code again, the limits are checked every time there is a call into CVodesIntegrator::integrate(double tout) (where tout is based on the estimated advance limit). It is true that CVODES takes several hundred steps internally, but I'm not sure how to check the advance limits from within vendored code.