Skip to content

Add mass matrix support for FIRK and MIRK#536

Draft
ChrisRackauckas-Claude wants to merge 1 commit into
SciML:masterfrom
ChrisRackauckas-Claude:pr521-mass-matrix-current
Draft

Add mass matrix support for FIRK and MIRK#536
ChrisRackauckas-Claude wants to merge 1 commit into
SciML:masterfrom
ChrisRackauckas-Claude:pr521-mass-matrix-current

Conversation

@ChrisRackauckas-Claude

Copy link
Copy Markdown
Contributor

This PR should be ignored until reviewed by @ChrisRackauckas.

Continues and rebases the mass matrix support from #521 onto current master without pushing to the contributor branch.

What changed:

  • Adds mass matrix handling for FIRK and MIRK collocation paths.
  • Adds shared Core helpers for mass-matrix stage application and algebraic constraint projection.
  • Rejects mesh adaptivity for DAE solves where the mass-matrix treatment is not supported.
  • Adds focused FIRK and MIRK DAE regression tests.

Local validation run on Julia +1:

  • timeout 3600 ~/.juliaup/bin/julia +1 --project=lib/BoundaryValueDiffEqMIRK -e 'import Pkg; Pkg.instantiate(); include("lib/BoundaryValueDiffEqMIRK/test/Core/dae_tests.jl")'
    • Simple index-1 DAE: 30 passed
    • Ascher & Spiteri example problem 1: 8 passed
    • Mesh adaptivity is not supported for DAEs: 2 passed
  • timeout 3600 ~/.juliaup/bin/julia +1 --project=lib/BoundaryValueDiffEqFIRK -e 'import Pkg; Pkg.instantiate(); include("lib/BoundaryValueDiffEqFIRK/test/expanded/dae_tests.jl")'
    • Simple index-1 DAE: 24 passed
    • Ascher & Spiteri example problem 1: 4 passed
    • Mesh adaptivity is not supported for DAEs: 2 passed
  • timeout 3600 ~/.juliaup/bin/julia +1 --project=lib/BoundaryValueDiffEqFIRK -e 'import Pkg; Pkg.instantiate(); include("lib/BoundaryValueDiffEqFIRK/test/nested/dae_tests.jl")'
    • Simple index-1 DAE: 18 passed
    • Mesh adaptivity is not supported for DAEs: 2 passed
  • timeout 3600 ~/.juliaup/bin/julia +1 --project=.runic_env -m Runic --check --diff . passed with no diff; .runic_env was removed after the check.

Add mass matrix handling for FIRK and MIRK collocation paths, including algebraic constraint application and DAE adaptivity guards.

Add focused FIRK and MIRK DAE tests for simple index-1 problems and Ascher-Spiteri coverage.

Co-Authored-By: Sebastian Micluța-Câmpeanu <sebastian.mc95@proton.me>

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
@github-actions

github-actions Bot commented Jul 8, 2026

Copy link
Copy Markdown
Contributor

Benchmark Results (Julia v1.11)

Time benchmarks
master d2f3541... master / d2f3541...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2() 0.576 ± 0.016 s 0.546 ± 0.0071 s 1.05 ± 0.033
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3() 12.3 ± 1.3 ms 12.2 ± 1.6 ms 1.01 ± 0.17
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4() 2.75 ± 0.44 ms 2.5 ± 0.19 ms 1.1 ± 0.19
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5() 3.25 ± 0.47 ms 3.05 ± 0.28 ms 1.07 ± 0.18
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6() 1.61 ± 0.32 ms 1.45 ± 0.14 ms 1.11 ± 0.24
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false) 1.71 ± 0.45 ms 1.6 ± 0.22 ms 1.07 ± 0.32
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true) 3.6 ± 1.1 ms 3.3 ± 0.56 ms 1.09 ± 0.38
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false) 0.0578 ± 0.017 s 0.0394 ± 0.0053 s 1.47 ± 0.49
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true) 0.0824 ± 0.014 s 0.0798 ± 0.0097 s 1.03 ± 0.22
Simple Pendulum/IIP/Shooting(Tsit5()) 0.191 ± 0.092 ms 0.166 ± 0.081 ms 1.15 ± 0.8
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2() 0.823 ± 0.023 s 0.836 ± 0.019 s 0.984 ± 0.036
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3() 17.1 ± 8.9 ms 19.4 ± 7.1 ms 0.882 ± 0.56
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4() 3.38 ± 0.65 ms 3.65 ± 0.94 ms 0.926 ± 0.3
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5() 4.2 ± 0.75 ms 4.76 ± 0.87 ms 0.882 ± 0.23
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6() 1.9 ± 0.4 ms 2.13 ± 0.5 ms 0.892 ± 0.28
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false) 4.01 ± 3.3 ms 4.1 ± 3.3 ms 0.979 ± 1.1
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true) 8.29 ± 7.4 ms 8.87 ± 8.3 ms 0.934 ± 1.2
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false) 0.111 ± 0.062 s 0.134 ± 0.057 s 0.831 ± 0.58
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true) 0.195 ± 0.057 s 0.156 ± 0.0093 s 1.25 ± 0.37
Simple Pendulum/OOP/Shooting(Tsit5()) 0.622 ± 0.49 ms 0.62 ± 0.44 ms 1 ± 1.1
time_to_load 7.83 ± 0.36 s 7.2 ± 0.31 s 1.09 ± 0.069
Memory benchmarks
master d2f3541... master / d2f3541...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2() 0.389 M allocs: 0.0443 GB 0.389 M allocs: 0.0443 GB 1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3() 0.0436 M allocs: 4.88 MB 0.0436 M allocs: 4.88 MB 1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4() 15.7 k allocs: 1.65 MB 15.7 k allocs: 1.65 MB 1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5() 22.1 k allocs: 2.03 MB 22.1 k allocs: 2.03 MB 1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6() 12.8 k allocs: 1.05 MB 12.8 k allocs: 1.05 MB 1
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false) 25.5 k allocs: 1.82 MB 25.5 k allocs: 1.82 MB 1
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true) 0.0491 M allocs: 3.39 MB 0.0491 M allocs: 3.39 MB 1
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false) 0.553 M allocs: 0.0535 GB 0.553 M allocs: 0.0535 GB 1
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true) 0.833 M allocs: 0.0778 GB 0.833 M allocs: 0.0778 GB 1
Simple Pendulum/IIP/Shooting(Tsit5()) 4.63 k allocs: 0.224 MB 4.63 k allocs: 0.224 MB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2() 0.89 M allocs: 0.984 GB 0.89 M allocs: 0.984 GB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3() 0.0932 M allocs: 24.8 MB 0.0932 M allocs: 24.8 MB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4() 0.0324 M allocs: 3.96 MB 0.0324 M allocs: 3.96 MB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5() 0.045 M allocs: 4.99 MB 0.045 M allocs: 4.99 MB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6() 25.3 k allocs: 2.17 MB 25.3 k allocs: 2.17 MB 1
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false) 0.142 M allocs: 10.2 MB 0.142 M allocs: 10.2 MB 1
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true) 0.266 M allocs: 18.7 MB 0.266 M allocs: 18.7 MB 1
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false) 2.52 M allocs: 0.279 GB 2.52 M allocs: 0.279 GB 1
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true) 3.82 M allocs: 0.404 GB 3.82 M allocs: 0.404 GB 1
Simple Pendulum/OOP/Shooting(Tsit5()) 0.0373 M allocs: 1.69 MB 0.0373 M allocs: 1.69 MB 1
time_to_load 0.161 k allocs: 11.4 kB 0.161 k allocs: 11.4 kB 1

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