Skip to content

WIP: add fixed-step RKV87-IIa (not yet working) #2693

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

WIP: add fixed-step RKV87-IIa (not yet working) #2693

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ PrecompileTools.@compile_workload begin
solver_list = nothing
end

export Vern6, Vern7, Vern8, Vern9
export AutoVern6, AutoVern7, AutoVern8, AutoVern9
export Vern6, Vern7, Vern8, Vern9, RKV87
export AutoVern6, AutoVern7, AutoVern8, AutoVern9, RKV87

end
6 changes: 6 additions & 0 deletions lib/OrdinaryDiffEqVerner/src/alg_utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,4 +12,10 @@ alg_stability_size(alg::Vern7) = 4.6400
alg_stability_size(alg::Vern8) = 5.8641
alg_stability_size(alg::Vern9) = 4.4762

# ── Verner RKV87-IIa traits ────────────────────────────────────────
alg_order(::RKV87) = 8 # primary solution order
embedded_order(::RKV87) = 7 # will matter once adaptivity is added
isfsal(::RKV87) = false # pair is NOT FSAL
has_dense_output(::RKV87)= false # we haven’t wired the interpolant yet
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

that's not a trait, also it has a fallback to linear.


SciMLBase.has_lazy_interpolation(alg::Union{Vern6, Vern7, Vern8, Vern9}) = true
24 changes: 24 additions & 0 deletions lib/OrdinaryDiffEqVerner/src/algorithms.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,27 @@ function Vern8(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true)
Vern8(stage_limiter!, step_limiter!, False(), lazy)
end

# ------------------------------------------------------------------
# Verner RKV87-IIa : 13-stage explicit 8 (7) Runge–Kutta pair
# Fixed-step implementation (no adaptivity or dense output yet)
# ------------------------------------------------------------------

@doc explicit_rk_docstring(
"Verner’s “efficient” 13-stage 8/7 Runge–Kutta method (RKV87-IIa).\n" *
"This variant is wired **for fixed time steps**. Use it via\n\n" *
"```julia\nsol = solve(prob, RKV87(), dt = h)\n```",
"RKV87")
Base.@kwdef struct RKV87{StageLimiter, StepLimiter, Thread} <:
OrdinaryDiffEqConstantAlgorithm
stage_limiter!::StageLimiter = trivial_limiter!
step_limiter!::StepLimiter = trivial_limiter!
thread::Thread = False()
end
TruncatedStacktraces.@truncate_stacktrace RKV87 3

# Register the algorithm with the dispatcher
add_algorithm(RKV87(), RKV87ConstantCache, Vern87IIaTableau)

@doc explicit_rk_docstring(
"Verner's “Most Efficient” 9/8 Runge-Kutta method. (lazy9th order interpolant).",
"Vern9",
Expand Down Expand Up @@ -113,6 +134,7 @@ function Vern9(stage_limiter!, step_limiter! = trivial_limiter!; lazy = true)
Vern9(stage_limiter!, step_limiter!, False(), lazy)
end


"""
Automatic switching algorithm that can switch between the (non-stiff) `Vern6()` and `stiff_alg`.

Expand Down Expand Up @@ -153,3 +175,5 @@ To gain access to stiff algorithms you might have to install additional librarie
such as `OrdinaryDiffEqRosenbrock`.
"""
AutoVern9(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern9(lazy = lazy), alg; kwargs...)


13 changes: 8 additions & 5 deletions lib/OrdinaryDiffEqVerner/src/verner_caches.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -263,9 +263,12 @@ struct Vern9ConstantCache <: OrdinaryDiffEqConstantCache
lazy::Bool
end

function alg_cache(alg::Vern9, u, rate_prototype, ::Type{uEltypeNoUnits},
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
dt, reltol, p, calck,
::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
Vern9ConstantCache(alg.lazy)
# ------------------------------------------------------------------
# Cache for Verner RKV87-IIa (13 stage values + 1 work vector)
# ------------------------------------------------------------------
@cache RKV87ConstantCache begin
k1::K ; k2::K ; k3::K ; k4::K ; k5::K ; k6::K ; k7::K
k8::K ; k9::K ; k10::K ; k11::K ; k12::K ; k13::K
tmp::K # work vector for low-storage kernel
end
Comment on lines +266 to 273
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're missing the alg_cache dispatches and the definition of the working cache.


12 changes: 12 additions & 0 deletions lib/OrdinaryDiffEqVerner/src/verner_rk_perform_step.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1013,6 +1013,18 @@ end
end
end

# ------------------------------------------------------------------
# RKV87-IIa wrappers
# ------------------------------------------------------------------

initialize!(integrator, cache::RKV87ConstantCache) = nothing

# The tableau we built is an ExplicitRKTableau with 13 stages
perform_step!(integrator,
cache::RKV87ConstantCache,
tab::ExplicitRKTableau{13}) =
verner_lowstorage_step!(integrator, cache, tab)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function doesn't exist? Is this outputted by an LLM?


function initialize!(integrator, cache::Vern9Cache)
@unpack k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15, k16 = cache
@unpack k = integrator
Expand Down
157 changes: 157 additions & 0 deletions lib/OrdinaryDiffEqVerner/src/verner_tableaus.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
using StaticArrays
module OrdinaryDiffEqVerner
## Vern6
struct Vern6ExtraStages{T, T2}
c10::T2
Expand Down Expand Up @@ -2306,6 +2308,83 @@ function Vern8Tableau(T, T2)
extra, interp)
end

# --------------------------------------------------------------------
# Verner RKV87-IIa : 13-stage explicit 8 (7) pair (efficient family)
# --------------------------------------------------------------------
function Vern87IIaTableau(::Type{T}, ::Type{T2}) where {T, T2}
c = @SVector T[
0.0,
0.092662,
0.131223036175,
0.196834554263,
0.427173,
0.485972,
0.161915,
0.985468,
0.96269773486,
0.99626,
0.997947,
1.0,
1.0
]

A = @SMatrix T[
0 0 0 0 0 0 0 0 0 0 0 0 0;
0.092662 0 0 0 0 0 0 0 0 0 0 0 0;
0.0383074654825 0.0929155706929 0 0 0 0 0 0 0 0 0 0 0;
0.0492086385658 0 0.147625915697 0 0 0 0 0 0 0 0 0 0;
0.27430760857 0 -0.93198872031 1.08485411174 0 0 0 0 0 0 0 0 0;
0.0646185297094 0 0 0.268762921337 0.152590548954 0 0 0 0 0 0 0 0;
0.0718915581977 0 0 0.312986472775 0.140861770313 0.0559558792414 0 0 0 0 0 0 0;
0.159677185234 0 0 0.693559173307 -0.316843836922 0.420825919995 -0.0847501736448 0 0 0 0 0 0;
0.0591917679674 0 0 0 -2.23377950868 0.390944379945 -0.936828383822 58.4572000999 0.00589430972373 0 0 0 0;
-6.68986189932 0 0 -81.4427100405 13.36778825698 -4.47077763842 80.2332139216 -0.0131363833621 0.0117437830322 0 0 0 0;
-6.7888419558 0 0 -82.6563985593 13.5997392187 -4.65073689329 82.1513831957 -0.0135453129954 0.0121304908312 -0.000312288856455 0 0 0;
0.0514137795362 0 0 0 0 0.318131086582 0.458946584938 -0.400240655913 0.0569851840034 -0.0010554370029 0.0000499398326837 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0
]

b = @SVector T[
0.0462554315971,
0.0,
0.0,
0.0,
0.0,
0.370666616552,
0.259044082455,
-679.984146817,
49.8916112904,
10271.2352221373,
-14782.1966063569,
5141.37795361606,
0.0
]

b̂ = @SVector T[ # embedded order-7 weights
0.0463850423437,
0.0,
0.0,
0.0,
0.0,
0.372576768158,
0.258568549512,
-147.495076759,
23.8436271264,
347.426416673,
0.0,
0.0,
-223.452497401
]

return ExplicitRKTableau(c, A, b, b̂)
end

# Fast path for CompiledFloats (optional; fine to add later)
function Vern87IIaTableau(::Type{<:CompiledFloats}, ::Type{<:CompiledFloats})
Vern87IIaTableau(Float64, Float64)
end


## Vern9
struct Vern9ExtraStages{T, T2}
c17::T2
Expand Down Expand Up @@ -3891,4 +3970,82 @@ end
a1610, a1611, a1612, a1613, b1, b8, b9, b10, b11, b12, b13, b14, b15,
btilde1, btilde8, btilde9, btilde10, btilde11, btilde12, btilde13,
btilde14, btilde15, btilde16)


end

# --------------------------------------------------------------------
# Verner RKV87-IIa : 13-stage explicit 8 (7) pair (efficient family)
# --------------------------------------------------------------------
function Vern87IIaTableau(::Type{T}, ::Type{T2}) where {T, T2}
c = @SVector T[
0.0,
0.092662,
0.131223036175,
0.196834554263,
0.427173,
0.485972,
0.161915,
0.985468,
0.96269773486,
0.99626,
0.997947,
1.0,
1.0
]

A = @SMatrix T[
0 0 0 0 0 0 0 0 0 0 0 0 0;
0.092662 0 0 0 0 0 0 0 0 0 0 0 0;
0.0383074654825 0.0929155706929 0 0 0 0 0 0 0 0 0 0 0;
0.0492086385658 0 0.147625915697 0 0 0 0 0 0 0 0 0 0;
0.27430760857 0 -0.93198872031 1.08485411174 0 0 0 0 0 0 0 0 0;
0.0646185297094 0 0 0.268762921337 0.152590548954 0 0 0 0 0 0 0 0;
0.0718915581977 0 0 0.312986472775 0.140861770313 0.0559558792414 0 0 0 0 0 0 0;
0.159677185234 0 0 0.693559173307 -0.316843836922 0.420825919995 -0.0847501736448 0 0 0 0 0 0;
0.0591917679674 0 0 0 -2.23377950868 0.390944379945 -0.936828383822 58.4572000999 0.00589430972373 0 0 0 0;
-6.68986189932 0 0 -81.4427100405 13.36778825698 -4.47077763842 80.2332139216 -0.0131363833621 0.0117437830322 0 0 0 0;
-6.7888419558 0 0 -82.6563985593 13.5997392187 -4.65073689329 82.1513831957 -0.0135453129954 0.0121304908312 -0.000312288856455 0 0 0;
0.0514137795362 0 0 0 0 0.318131086582 0.458946584938 -0.400240655913 0.0569851840034 -0.0010554370029 0.0000499398326837 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0
]

b = @SVector T[
0.0462554315971,
0.0,
0.0,
0.0,
0.0,
0.370666616552,
0.259044082455,
-679.984146817,
49.8916112904,
10271.2352221373,
-14782.1966063569,
5141.37795361606,
0.0
]

b̂ = @SVector T[ # embedded order-7 weights
0.0463850423437,
0.0,
0.0,
0.0,
0.0,
0.372576768158,
0.258568549512,
-147.495076759,
23.8436271264,
347.426416673,
0.0,
0.0,
-223.452497401
]

return ExplicitRKTableau(c, A, b, b̂)
end

# Fast path for CompiledFloats (optional; fine to add later)
function Vern87IIaTableau(::Type{<:CompiledFloats}, ::Type{<:CompiledFloats})
Vern87IIaTableau(Float64, Float64)
end