advance limit handling to improve performance (#1453)

[wandadars]

This was just an attempt to try and resolve the issue discussed in #1453. This may be too invasive to be useful to the code base, but just wanted to open a pull request just in case.

Event-based limiter: Enforce advance limits with a CVODE root function so integration stops exactly when any component reaches its configured limit, preventing large overshoots between output points. Hooked via CVodeRootInit and handled CV_ROOT_RETURN.

Closes #1453

First example case from #1453

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()
# Show plot on screen
plt.show()

The plot from that case using this updated approach is:

I also ran the case from the original issue report with the unit test parameters.

Click to view code example 1

import cantera as ct
import matplotlib.pyplot as plt


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)

    # Unit test parameters
    tEnd = 0.01
    tStep = 1e-3
    nSteps = 0

    #tEnd = 0.01
    #tStep = 1.1e-4
    #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)

fig, axes = plt.subplots(
    nrows=4, ncols=1, figsize=(6, 10), sharex=True,
    gridspec_kw={'height_ratios': [2, 1, 1, 1], 'hspace': 0.25}
)
ax_main, ax_fine, ax_coarse, ax_baseline = axes

# Main plot: fine and coarse together
ax_main.plot(states_fine.t, states_fine('H2').Y, '-', color='orange', label='fine (line)')
ax_main.plot(states_fine.t, states_fine('H2').Y, '.', color='orange', label='fine (points)')
ax_main.plot(states_coarse.t, states_coarse('H2').Y, '--', color='blue', label='coarse (dashed)')
ax_main.plot(states_coarse.t, states_coarse('H2').Y, 'o', color='blue', label='coarse (points)')
ax_main.plot(states_baseline.t, states_baseline('H2').Y, '-', color='green', label='baseline')
ax_main.plot(states_baseline.t, states_baseline('H2').Y, '.', color='green')
ax_main.set_ylabel('Y(H2)')
ax_main.legend(loc='best')

# Subplot 1: fine only
ax_fine.plot(states_fine.t, states_fine('H2').Y, '-', color='orange', label='fine')
ax_fine.plot(states_fine.t, states_fine('H2').Y, '.', color='orange')
ax_fine.set_ylabel('Y(H2)')
ax_fine.legend(loc='best')

# Subplot 2: coarse only
ax_coarse.plot(states_coarse.t, states_coarse('H2').Y, '--', color='blue', label='coarse')
ax_coarse.plot(states_coarse.t, states_coarse('H2').Y, 'o', color='blue')
ax_coarse.set_ylabel('Y(H2)')
ax_coarse.legend(loc='best')

# Subplot 3: baseline only
ax_baseline.plot(states_baseline.t, states_baseline('H2').Y, '-', color='green', label='baseline')
ax_baseline.plot(states_baseline.t, states_baseline('H2').Y, '.', color='green')
ax_baseline.set_ylabel('Y(H2)')
ax_baseline.set_xlabel('time [s]')
ax_baseline.legend(loc='best')

plt.show()

This was using the default values from the unit test

And this was using the tighter tolerance values from the issue discussion.

case: limit=None, dy_max=0.00551594315108739, nSteps=90
case: limit=0.01, dy_max=0.00551594315108739, nSteps=90
case: limit=0.001, dy_max=0.0010000000000554895, nSteps=95

For reference, this is the behavior of the Cantera version (3.1.0) running the same script with the tighter tolerances.

case: limit=None, dy_max=0.0055159431511252995, nSteps=90
case: limit=0.01, dy_max=0.005365298679090868, nSteps=93
case: limit=0.001, dy_max=0.005365298679090868, nSteps=93

[speth]

Thanks for looking into this, @wandadars. I think using the built-in CVODES functionality for this is a good approach to consider. I just wanted to share a couple of suggestions that might help with the implementation.

[speth]

I suspect you can avoid breaking the abstraction between CVodesIntegrator and ReactorNet by making the call to CVodeRootInit optional, i.e. triggered only if a call to advance with limits applied is made. The CVodes documentation seems to suggest that repeated calls to CVodeRootInit changing the number of root functions to evaluate (including setting it to zero) are allowed.

[speth]

I think the current test failure (integrator hanging) stems from this. There is a tricky issue with both CVodesIntegrator and ReactorNet having a m_time member, where the former is shadowed by the latter. What I think is happening here is that only the former is being correctly updated to the user-specified time with the interpolated solution while the latter ends up having the time reached by the integrator internally on the last step.