Skip to content

Add mass matrix support for FIRK and MIRK#521

Open
SebastianM-C wants to merge 2 commits into
SciML:masterfrom
SebastianM-C:smc/mass_matrix
Open

Add mass matrix support for FIRK and MIRK#521
SebastianM-C wants to merge 2 commits into
SciML:masterfrom
SebastianM-C:smc/mass_matrix

Conversation

@SebastianM-C

@SebastianM-C SebastianM-C commented Jul 2, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

via unprojected collocation for index-1 DAEs based on Ascher & Spiteri (1994)

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

Add any other context about the problem here.

@SebastianM-C @claude
Add mass matrix support for FIRK and MIRK
94a7cb6 
via unprojected collocation for index-1 DAEs based on Ascher & Spiteri (1994)

Co-authored-by: Claude Opus 4.6 <noreply@anthropic.com>
@github-actions

github-actions Bot commented Jul 2, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor

Benchmark Results (Julia v1.11)

Time benchmarks
master bd2fc6a... master / bd2fc6a...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2() 0.658 ± 0.0082 s 0.653 ± 0.013 s 1.01 ± 0.024
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3() 13.2 ± 0.25 ms 13.2 ± 0.31 ms 1 ± 0.031
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4() 2.81 ± 0.073 ms 2.8 ± 0.078 ms 1 ± 0.038
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5() 3.45 ± 0.093 ms 3.42 ± 0.11 ms 1.01 ± 0.042
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6() 1.65 ± 0.05 ms 1.65 ± 0.055 ms 0.998 ± 0.045
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false) 1.63 ± 0.098 ms 1.61 ± 0.096 ms 1.01 ± 0.086
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true) 3.2 ± 0.13 ms 3.18 ± 0.15 ms 1.01 ± 0.061
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false) 0.0405 ± 0.0035 s 0.0409 ± 0.0037 s 0.992 ± 0.12
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true) 0.0651 ± 0.002 s 0.0646 ± 0.00097 s 1.01 ± 0.034
Simple Pendulum/IIP/Shooting(Tsit5()) 0.178 ± 0.083 ms 0.177 ± 0.083 ms 1 ± 0.67
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2() 0.808 ± 0.0079 s 0.797 ± 0.014 s 1.01 ± 0.02
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3() 21.9 ± 5.9 ms 17.1 ± 6.2 ms 1.28 ± 0.58
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4() 3.42 ± 0.13 ms 3.43 ± 0.21 ms 0.997 ± 0.073
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5() 4.15 ± 0.2 ms 4.14 ± 0.21 ms 1 ± 0.07
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6() 1.96 ± 0.12 ms 1.94 ± 0.13 ms 1.01 ± 0.091
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false) 3.63 ± 0.9 ms 3.61 ± 1.8 ms 1.01 ± 0.57
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true) 6.93 ± 6.5 ms 6.84 ± 6.7 ms 1.01 ± 1.4
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false) 0.0977 ± 0.0034 s 0.0971 ± 0.0034 s 1.01 ± 0.05
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true) 0.151 ± 0.0055 s 0.153 ± 0.007 s 0.988 ± 0.058
Simple Pendulum/OOP/Shooting(Tsit5()) 0.643 ± 0.051 ms 0.633 ± 0.045 ms 1.02 ± 0.11
time_to_load 6.56 ± 0.01 s 6.54 ± 0.024 s 1 ± 0.004
Memory benchmarks
master bd2fc6a... master / bd2fc6a...
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.87 MB 0.0436 M allocs: 4.87 MB 1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4() 15.7 k allocs: 1.64 MB 15.7 k allocs: 1.64 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.04 MB 12.8 k allocs: 1.04 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.95 MB 0.0324 M allocs: 3.95 MB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5() 0.045 M allocs: 4.98 MB 0.045 M allocs: 4.98 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.159 k allocs: 11.2 kB 0.159 k allocs: 11.2 kB 1

@SebastianM-C @claude
Error on adaptive solves for DAEs and add index-1 DAE tests
bd2fc6a 
Defect-based mesh adaptivity cannot converge for DAEs since the
collocation interpolant is inaccurate for algebraic variables, so
requesting it now throws an ArgumentError pointing users to
`adaptive = false` or BoundaryValueDiffEqAscher.jl.

Also fix the nested FIRK out-of-place collocation kernels calling the
in-place algebraic constraint helper on the user's out-of-place f.

Tests cover MIRK2-6, RadauIIa and LobattoIIIc (expanded and nested) on
a simple index-1 DAE and on example problem 1 from Ascher & Spiteri
(1994), plus the adaptivity error. LobattoIIIa/IIIb are excluded: their
tableau structure cannot determine the algebraic stage components.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@SebastianM-C SebastianM-C marked this pull request as ready for review July 3, 2026 14:11
@ChrisRackauckas-Claude ChrisRackauckas-Claude mentioned this pull request Jul 8, 2026
Draft
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

1 participant

@SebastianM-C