Skip to content

Conversation

@ChrisRackauckas-Claude
Copy link
Contributor

Summary

  • Fixes issue Stabilized RK methods do not work with Unitful #2858 where stabilized RK methods (ESERK4, ESERK5, SERK2, ROCK2, ROCK4, RKC) failed when time was given as a Unitful quantity
  • Root cause: The maxeig! power iteration function in rkc_utils.jl mixed quantities with rate units [u]/[t] and state units [u], causing dimensional incompatibility
  • Key fixes:
    • Track whether z has rate units and multiply by dt when converting to state units
    • Use DiffEqBase.value() to strip units from dt and dtmax when computing unitless eigenvalue estimates
    • Fix quotient formula in power iteration loop to maintain constant perturbation size
    • Fix sqrt(abs(dt) * integrator.eigen_est) in perform_step to use unitless dt value

Test plan

  • All 68 tests in OrdinaryDiffEqStabilizedRK pass locally
  • Verified Unitful compatibility with exact test case from issue Stabilized RK methods do not work with Unitful #2858:
    using OrdinaryDiffEq, Unitful
    function lorenz!(du, u, p, t)
        du[1] = 10.0*(u[2]-u[1])
        du[2] = u[1]*(28.0-u[3]) - u[2]
        du[3] = u[1]*u[2] - (8/3)*u[3]
    end
    u0 = [1.0; 0.0; 0.0]
    tspan = (0.0u"s", 100.0u"s")
    prob = ODEProblem(lorenz!, u0, tspan)
    sol = solve(prob, ROCK4())  # Previously threw DimensionError
  • CI tests

🤖 Generated with Claude Code

Co-Authored-By: Claude Opus 4.5 [email protected]

This commit fixes issue SciML#2858 where stabilized RK methods (ESERK4, ESERK5,
SERK2, ROCK2, ROCK4, RKC) failed when time was given as a Unitful quantity.

The root cause was dimensional incompatibility in the maxeig! power iteration
function and various places in the perform_step implementations:

1. In maxeig!, fsalfirst has rate units [u]/[t] while z needs state units [u]
   for the perturbation z - uprev to be valid. Fixed by multiplying by dt when
   converting rate to state, and adjusting the normalization formula to maintain
   correct perturbation size.

2. In perform_step for mdeg calculation, sqrt(abs(dt) * eigen_est) mixed units.
   Fixed by using DiffEqBase.value(abs(dt)) to get unitless dt.

3. In ROCK2/ROCK4, internalnorm was incorrectly applied to time values.
   Fixed by using time differences directly.

4. In RKC error estimate, the formula mixed rate and state units.
   Fixed by using the last derivative evaluation (feval) instead of gprev.

All existing tests pass and Unitful compatibility is verified for both
in-place and out-of-place problem formulations.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Opus 4.5 <[email protected]>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants