advance limit handling to improve performance (#1453) by wandadars · Pull Request #1976 · Cantera/cantera

wandadars · 2025-09-12T19:38:57Z

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.

src/numerics/CVodesIntegrator.cpp

src/zeroD/ReactorNet.cpp

wandadars · 2025-11-01T02:52:33Z

Thanks @speth . I totally spaced out and forgot that I had opened this PR. I'll look at this again this weekend and see if I can clean it up.

codecov · 2025-11-02T04:34:20Z

Codecov Report

❌ Patch coverage is 70.64220% with 32 lines in your changes missing coverage. Please review.
✅ Project coverage is 76.38%. Comparing base (61d5c5b) to head (b29bba6).
⚠️ Report is 16 commits behind head on main.

Files with missing lines	Patch %	Lines
src/numerics/FuncEval.cpp	14.28%	18 Missing ⚠️
src/zeroD/ReactorNet.cpp	85.96%	2 Missing and 6 partials ⚠️
src/numerics/CVodesIntegrator.cpp	82.60%	2 Missing and 2 partials ⚠️
include/cantera/numerics/FuncEval.h	66.66%	1 Missing ⚠️
include/cantera/numerics/Integrator.h	0.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1976      +/-   ##
==========================================
- Coverage   76.41%   76.38%   -0.03%     
==========================================
  Files         463      463              
  Lines       54952    55032      +80     
  Branches     9308     9317       +9     
==========================================
+ Hits        41990    42038      +48     
- Misses       9831     9863      +32     
  Partials     3131     3131

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

ischoegl

Thanks! It's great to see that there's actually a way to handle this. My only request would be to add a test so this is covered.

I hope to see this in the final version of 3.2 ...

speth

Thanks, @wandadars, this looks really good. I was worried that this would introduce a significant performance penalty due to the iteration required to exactly hit the target advance limit, but I think what actually happens is that this iteration is done entirely on the interpolated solution vector and as such doesn't require any additional evaluations of the ODE RHS function.

Addendum: I also agree with @ischoegl's suggestion on improving the testing of this. Probably the biggest gap in the coverage is on the verbose logging side. You can test this in pytest using the capsys fixture; there are a couple of instances in the test suite now that use this that you can use as examples.

include/cantera/numerics/FuncEval.h

include/cantera/numerics/Integrator.h

include/cantera/zeroD/ReactorNet.h

src/zeroD/ReactorNet.cpp

src/numerics/CVodesIntegrator.cpp

src/zeroD/ReactorNet.cpp

wandadars · 2025-11-13T06:27:45Z

Thanks for the comments @speth and @ischoegl . I believe all the comments have been addressed and tests have been added/updated regarding the advance limits.

speth

Thanks for the updates, @wandadars. This looks good to me.

I added the "no throw" wrapper around evalRootFunctions that I commented on but may not have been completely clear about and updated a few comments.

@ischoegl, since I've now modified this PR, will you please review/approve?

ischoegl

Thank you, @wandadars (and @speth). As the original author of the advance limits capability, I am very happy to see this being made more functional.

ischoegl added the Reactors label Sep 13, 2025

speth reviewed Sep 15, 2025

View reviewed changes

src/numerics/CVodesIntegrator.cpp Outdated Show resolved Hide resolved

src/zeroD/ReactorNet.cpp Show resolved Hide resolved

wandadars force-pushed the advance_limits_improvements branch 2 times, most recently from d476b5c to 65704b2 Compare November 2, 2025 03:49

ischoegl reviewed Nov 4, 2025

View reviewed changes

speth requested changes Nov 5, 2025

View reviewed changes

wandadars and others added 7 commits November 13, 2025 03:26

updates to advance limit handling to improve performance

cf6a9d2

addressing review comments

db85f0f

updated comments for new functions

bae8bc5

fixed comment formatting

37e1784

fixed issues with Sundials callbacks

bd5e134

addressing review comments

c39f912

added/updated tests involving advance limits

57028f5

wandadars force-pushed the advance_limits_improvements branch from 02b8fee to 57028f5 Compare November 13, 2025 08:26

speth added 2 commits November 14, 2025 11:00

[Numerics] Handle possible exceptions in evalRootFunctions

aa67d77

[Doc] Minor adjustments to advance limit docs

b29bba6

speth approved these changes Nov 14, 2025

View reviewed changes

speth requested a review from ischoegl November 14, 2025 16:16

ischoegl approved these changes Nov 14, 2025

View reviewed changes

ischoegl merged commit d1d8536 into Cantera:main Nov 14, 2025
58 of 61 checks passed

Uh oh!

Conversation

wandadars commented Sep 12, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

speth left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

wandadars commented Nov 1, 2025

Uh oh!

codecov bot commented Nov 2, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

ischoegl left a comment

Choose a reason for hiding this comment

Uh oh!

speth left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

wandadars commented Nov 13, 2025

Uh oh!

speth left a comment

Choose a reason for hiding this comment

Uh oh!

ischoegl left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

wandadars commented Sep 12, 2025 •

edited

Loading

codecov bot commented Nov 2, 2025 •

edited

Loading

speth left a comment •

edited

Loading