From 937eb1a9952be34d80d9694d7274fa622b2f1ff9 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 14 Feb 2025 14:36:04 -0500 Subject: [PATCH 01/32] init --- lib/ImplicitDiscreteSolve/Project.toml | 4 +++ .../src/ImplicitDiscreteSolve.jl | 31 +++++++++++++++++++ lib/ImplicitDiscreteSolve/src/solve.jl | 6 ++++ 3 files changed, 41 insertions(+) create mode 100644 lib/ImplicitDiscreteSolve/Project.toml create mode 100644 lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl create mode 100644 lib/ImplicitDiscreteSolve/src/solve.jl diff --git a/lib/ImplicitDiscreteSolve/Project.toml b/lib/ImplicitDiscreteSolve/Project.toml new file mode 100644 index 000000000..edee1e00a --- /dev/null +++ b/lib/ImplicitDiscreteSolve/Project.toml @@ -0,0 +1,4 @@ +name = "ImplicitDiscreteSolve" +uuid = "3263718b-31ed-49cf-8a0f-35a466e8af96" +authors = ["vyudu "] +version = "0.1.0" diff --git a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl new file mode 100644 index 000000000..a6db0f936 --- /dev/null +++ b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl @@ -0,0 +1,31 @@ +module ImplicitDiscreteSolve + +using SciMLBase: AbstractNonlinearProblem +using SciMLBase +using NonlinearSolveBase +using SymbolicIndexingInterface +using LinearAlgebra +using ADTypes +using TaylorDiff +using DocStringExtensions +import CommonSolve +import DifferentiationInterface as DI + +using ConcreteStructs: @concrete + +""" + IteratedNonlinearSolve(; nlsolvealg, autodiff = true, kwargs...) + +This algorithm is a solver for ImplicitDiscreteSystems. +""" +@concrete struct IteratedNonlinearSolve <: NonlinearSolveBase.AbstractNonlinearSolveAlgorithm + nlsolvealg + autodiff + kwargs +end + +export IteratedNonlinearSolve + +include("solve.jl") + +end diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl new file mode 100644 index 000000000..589a78151 --- /dev/null +++ b/lib/ImplicitDiscreteSolve/src/solve.jl @@ -0,0 +1,6 @@ + +# make a nonlinear problem and solve at every timestep. + +function CommonSolve.solve(prob::ImplicitDiscreteProblem) + +end From 55cecfe2518fe12eaf745e97a217f2747a82e550 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 19 Feb 2025 17:33:47 -0800 Subject: [PATCH 02/32] begin implementing NLSolve --- lib/ImplicitDiscreteSolve/Project.toml | 34 ++++++++++++++ .../src/ImplicitDiscreteSolve.jl | 22 ++++----- lib/ImplicitDiscreteSolve/src/cache.jl | 35 ++++++++++++++ lib/ImplicitDiscreteSolve/src/solve.jl | 46 +++++++++++++++++-- lib/ImplicitDiscreteSolve/test/runtests.jl | 26 +++++++++++ 5 files changed, 149 insertions(+), 14 deletions(-) create mode 100644 lib/ImplicitDiscreteSolve/src/cache.jl create mode 100644 lib/ImplicitDiscreteSolve/test/runtests.jl diff --git a/lib/ImplicitDiscreteSolve/Project.toml b/lib/ImplicitDiscreteSolve/Project.toml index edee1e00a..527a2e957 100644 --- a/lib/ImplicitDiscreteSolve/Project.toml +++ b/lib/ImplicitDiscreteSolve/Project.toml @@ -2,3 +2,37 @@ name = "ImplicitDiscreteSolve" uuid = "3263718b-31ed-49cf-8a0f-35a466e8af96" authors = ["vyudu "] version = "0.1.0" + +[deps] +ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" +CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" +DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0" +OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" +OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" +SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" +UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" + +[compat] +ADTypes = "1.13.0" +CommonSolve = "0.2.4" +DiffEqBase = "6.164.1" +DifferentiationInterface = "0.6.42" +LinearAlgebra = "1.11.0" +NonlinearSolveBase = "1.5.0" +OrdinaryDiffEqCore = "1.18.1" +OrdinaryDiffEqSDIRK = "1.2.0" +Reexport = "1.2.2" +SciMLBase = "2.74.1" +SimpleNonlinearSolve = "2.1.0" +SymbolicIndexingInterface = "0.3.37" +UnPack = "1.0.2" + +[extras] +OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" +SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" diff --git a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl index a6db0f936..4ffcf234b 100644 --- a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl +++ b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl @@ -6,26 +6,26 @@ using NonlinearSolveBase using SymbolicIndexingInterface using LinearAlgebra using ADTypes -using TaylorDiff -using DocStringExtensions +using UnPack +import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, initialize!, perform_step! import CommonSolve import DifferentiationInterface as DI -using ConcreteStructs: @concrete +using Reexport +@reexport using DiffEqBase """ - IteratedNonlinearSolve(; nlsolvealg, autodiff = true, kwargs...) + IDSolve(alg; autodiff = true, kwargs...) -This algorithm is a solver for ImplicitDiscreteSystems. +Solver for `ImplicitDiscreteSystems`. `alg` is the NonlinearSolve algorithm that is used to solve for the next timestep at each step. """ -@concrete struct IteratedNonlinearSolve <: NonlinearSolveBase.AbstractNonlinearSolveAlgorithm - nlsolvealg - autodiff - kwargs +struct IDSolve{algType} <: OrdinaryDiffEqAlgorithm + nlsolve::algType end -export IteratedNonlinearSolve - +include("cache.jl") include("solve.jl") +export IDSolve + end diff --git a/lib/ImplicitDiscreteSolve/src/cache.jl b/lib/ImplicitDiscreteSolve/src/cache.jl new file mode 100644 index 000000000..4eef79141 --- /dev/null +++ b/lib/ImplicitDiscreteSolve/src/cache.jl @@ -0,0 +1,35 @@ +mutable struct ImplicitDiscreteState{uType, pType, tType} + u::Vector{uType} + p::Union{Nothing, Vector{pType}} + t_next::tType +end + +mutable struct IDSolveCache{uType} <: OrdinaryDiffEqMutableCache + u::uType + uprev::uType + state::ImplicitDiscreteState + prob::Union{Nothing, AbstractNonlinearProblem} +end + +function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + + state = ImplicitDiscreteState(similar(u), similar(p), t) + IDSolveCache(u, uprev, state, nothing) +end + +isdiscretecache(cache::IDSolveCache) = true + +function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + + state = ImplicitDiscreteState(similar(u), similar(p), t) + IDSolveCache(u, uprev, state, nothing) +end + +isfsal(alg::IDSolve) = false +get_fsalfirstlast(cache::IDSolveCache, rate_prototype) = (nothing, nothing) diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl index 589a78151..9e409bf0d 100644 --- a/lib/ImplicitDiscreteSolve/src/solve.jl +++ b/lib/ImplicitDiscreteSolve/src/solve.jl @@ -1,6 +1,46 @@ +# Remake the nonlinear problem, then update +function perform_step!(integrator, cache::IDSolveCache, repeat_step = false) + @unpack alg, u, uprev, dt, t, f, p = integrator + nlsolve = alg.nlsolve + @unpack state, prob = cache + state.u .= uprev + state.t_next = t + @show state + prob = remake(prob, p = state) -# make a nonlinear problem and solve at every timestep. + u = solve(prob, nlsolve) + any(isnan, u) && (integrator.sol.retcode = SciMLBase.ReturnCode.Failure) + integrator.u = u +end + +function initialize!(integrator, cache::IDSolveCache) + cache.state.u .= integrator.u + cache.state.p .= integrator.p + cache.state.t_next = integrator.t + f = integrator.f -function CommonSolve.solve(prob::ImplicitDiscreteProblem) - + _f = if isinplace(f) + (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t_next) + else + (u_next, p) -> f(u_next, p.u, p.p, p.t_next) + end + + prob = if isinplace(f) + NonlinearProblem{true}(_f, cache.state.u, cache.state) + else + NonlinearProblem{false}(_f, cache.state.u, cache.state) + end + cache.prob = prob end + +#### unnecessary +#function DiffEqBase.__init(prob::ImplicitDiscreteProblem, alg) +# f = prob.f +# t_i = prob.tspan[1] +# u0 = state_values(prob) +# p = parameter_values(prob) +# +# _f(resid, u_next, p) = f(resid, u_next, p.u, p.p, p.t) +# state = ImplicitDiscreteState(u0, p, t_i) +# nlprob = NonlinearProblem(_f, u0, state) +#end diff --git a/lib/ImplicitDiscreteSolve/test/runtests.jl b/lib/ImplicitDiscreteSolve/test/runtests.jl new file mode 100644 index 000000000..e9c4f3069 --- /dev/null +++ b/lib/ImplicitDiscreteSolve/test/runtests.jl @@ -0,0 +1,26 @@ +#runtests +using ImplicitDiscreteSolve +using OrdinaryDiffEqCore +using OrdinaryDiffEqSDIRK +using SimpleNonlinearSolve + +# Test implicit Euler using ImplicitDiscreteProblem +@testset "Implicit Discrete System" begin + function lotkavolterra(u, p, t) + [1.5*u[1] - u[1]*u[2], -3.0*u[2] + u[1]*u[2]] + end + + function f!(resid, u_next, u, p, t) + @. resid = u_next - 0.01*lotkavolterra(u_next, p, t) + end + u0 = [1., 1.] + tspan = (0., 1.) + + idprob = ImplicitDiscreteProblem(f!, u0, tspan, []; dt = 0.01) + idsol = solve(idprob, IDSolve(SimpleNewtonRaphson())) + + oprob = ODEProblem(lotkavolterra, u0, tspan) + osol = solve(oprob, ImplicitEuler()) + + @test idsol[end] ≈ osol[end] +end From f54da743eb9232c9a958c369b9063a21398cc153 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 20 Feb 2025 10:01:46 -0800 Subject: [PATCH 03/32] basic implementation working --- .../src/ImplicitDiscreteSolve.jl | 3 ++- lib/ImplicitDiscreteSolve/src/alg_utils.jl | 19 +++++++++++++++++++ lib/ImplicitDiscreteSolve/src/solve.jl | 13 ------------- lib/ImplicitDiscreteSolve/test/runtests.jl | 11 +++++++---- 4 files changed, 28 insertions(+), 18 deletions(-) create mode 100644 lib/ImplicitDiscreteSolve/src/alg_utils.jl diff --git a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl index 4ffcf234b..8fa8550d8 100644 --- a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl +++ b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl @@ -7,7 +7,7 @@ using SymbolicIndexingInterface using LinearAlgebra using ADTypes using UnPack -import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, initialize!, perform_step! +import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal, initialize!, perform_step!, isdiscretecache, isdiscretealg, alg_order, beta2_default, beta1_default, dt_required import CommonSolve import DifferentiationInterface as DI @@ -25,6 +25,7 @@ end include("cache.jl") include("solve.jl") +include("alg_utils.jl") export IDSolve diff --git a/lib/ImplicitDiscreteSolve/src/alg_utils.jl b/lib/ImplicitDiscreteSolve/src/alg_utils.jl new file mode 100644 index 000000000..5c25c2c21 --- /dev/null +++ b/lib/ImplicitDiscreteSolve/src/alg_utils.jl @@ -0,0 +1,19 @@ +function SciMLBase.isautodifferentiable(alg::IDSolve) + true +end +function SciMLBase.allows_arbitrary_number_types(alg::IDSolve) + true +end +function SciMLBase.allowscomplex(alg::IDSolve) + true +end + +SciMLBase.isdiscrete(alg::IDSolve) = true + +isfsal(alg::IDSolve) = false +alg_order(alg::IDSolve) = 0 +beta2_default(alg::IDSolve) = 0 +beta1_default(alg::IDSolve, beta2) = 0 + +dt_required(alg::IDSolve) = false +isdiscretealg(alg::IDSolve) = true diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl index 9e409bf0d..cf76a745e 100644 --- a/lib/ImplicitDiscreteSolve/src/solve.jl +++ b/lib/ImplicitDiscreteSolve/src/solve.jl @@ -5,7 +5,6 @@ function perform_step!(integrator, cache::IDSolveCache, repeat_step = false) @unpack state, prob = cache state.u .= uprev state.t_next = t - @show state prob = remake(prob, p = state) u = solve(prob, nlsolve) @@ -32,15 +31,3 @@ function initialize!(integrator, cache::IDSolveCache) end cache.prob = prob end - -#### unnecessary -#function DiffEqBase.__init(prob::ImplicitDiscreteProblem, alg) -# f = prob.f -# t_i = prob.tspan[1] -# u0 = state_values(prob) -# p = parameter_values(prob) -# -# _f(resid, u_next, p) = f(resid, u_next, p.u, p.p, p.t) -# state = ImplicitDiscreteState(u0, p, t_i) -# nlprob = NonlinearProblem(_f, u0, state) -#end diff --git a/lib/ImplicitDiscreteSolve/test/runtests.jl b/lib/ImplicitDiscreteSolve/test/runtests.jl index e9c4f3069..e8bb2b601 100644 --- a/lib/ImplicitDiscreteSolve/test/runtests.jl +++ b/lib/ImplicitDiscreteSolve/test/runtests.jl @@ -6,15 +6,18 @@ using SimpleNonlinearSolve # Test implicit Euler using ImplicitDiscreteProblem @testset "Implicit Discrete System" begin - function lotkavolterra(u, p, t) + function lotkavolterra(u, p, t) [1.5*u[1] - u[1]*u[2], -3.0*u[2] + u[1]*u[2]] end function f!(resid, u_next, u, p, t) - @. resid = u_next - 0.01*lotkavolterra(u_next, p, t) + lv = lotkavolterra(u_next, p, t) + resid[1] = u_next[1] - u[1] - 0.01*lv[1] + resid[2] = u_next[2] - u[2] - 0.01*lv[2] + nothing end u0 = [1., 1.] - tspan = (0., 1.) + tspan = (0., 0.5) idprob = ImplicitDiscreteProblem(f!, u0, tspan, []; dt = 0.01) idsol = solve(idprob, IDSolve(SimpleNewtonRaphson())) @@ -22,5 +25,5 @@ using SimpleNonlinearSolve oprob = ODEProblem(lotkavolterra, u0, tspan) osol = solve(oprob, ImplicitEuler()) - @test idsol[end] ≈ osol[end] + @test isapprox(idsol[end], osol[end], atol = 0.01) end From ea68831a9b5e67cbc599aef2554a6366f929bbe4 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 20 Feb 2025 10:10:55 -0800 Subject: [PATCH 04/32] bump Project.toml --- lib/ImplicitDiscreteSolve/Project.toml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/lib/ImplicitDiscreteSolve/Project.toml b/lib/ImplicitDiscreteSolve/Project.toml index 527a2e957..9e1096534 100644 --- a/lib/ImplicitDiscreteSolve/Project.toml +++ b/lib/ImplicitDiscreteSolve/Project.toml @@ -11,10 +11,8 @@ DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" -OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" @@ -36,3 +34,6 @@ UnPack = "1.0.2" [extras] OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" + +[targets] +test = ["OrdinaryDiffEqSDIRK", "SimpleNonlinearSolve"] From 1e14e64e44ca2bdb0e5bc4acb01d3e1ff872c69a Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 26 Feb 2025 10:15:59 -0500 Subject: [PATCH 05/32] refactor: refactor to --- .../src/SimpleImplicitDiscreteSolve.jl | 31 +++++++++++++++ .../src/alg_utils.jl | 19 +++++++++ lib/SimpleImplicitDiscreteSolve/src/cache.jl | 39 +++++++++++++++++++ lib/SimpleImplicitDiscreteSolve/src/solve.jl | 32 +++++++++++++++ 4 files changed, 121 insertions(+) create mode 100644 lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl create mode 100644 lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl create mode 100644 lib/SimpleImplicitDiscreteSolve/src/cache.jl create mode 100644 lib/SimpleImplicitDiscreteSolve/src/solve.jl diff --git a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl new file mode 100644 index 000000000..29f61cf66 --- /dev/null +++ b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl @@ -0,0 +1,31 @@ +module ImplicitDiscreteSolve + +using SciMLBase: AbstractNonlinearProblem +using SciMLBase +using NonlinearSolveBase +using SimpleNonlinearSolve +using SymbolicIndexingInterface +using LinearAlgebra +using ADTypes +using UnPack +import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal, initialize!, perform_step!, isdiscretecache, isdiscretealg, alg_order, beta2_default, beta1_default, dt_required +import CommonSolve +import DifferentiationInterface as DI + +using Reexport +@reexport using DiffEqBase + +""" + IDSolve(alg; autodiff = true, kwargs...) + +Solver for `ImplicitDiscreteSystems`. `alg` is the NonlinearSolve algorithm that is used to solve for the next timestep at each step. +""" +struct SimpleIDSolve{algType} <: OrdinaryDiffEqAlgorithm end + +include("cache.jl") +include("solve.jl") +include("alg_utils.jl") + +export SimpleIDSolve + +end diff --git a/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl b/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl new file mode 100644 index 000000000..b4eff1bc1 --- /dev/null +++ b/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl @@ -0,0 +1,19 @@ +function SciMLBase.isautodifferentiable(alg::SimpleIDSolve) + true +end +function SciMLBase.allows_arbitrary_number_types(alg::SimpleIDSolve) + true +end +function SciMLBase.allowscomplex(alg::SimpleIDSolve) + true +end + +SciMLBase.isdiscrete(alg::SimpleIDSolve) = true + +isfsal(alg::SimpleIDSolve) = false +alg_order(alg::SimpleIDSolve) = 0 +beta2_default(alg::SimpleIDSolve) = 0 +beta1_default(alg::SimpleIDSolve, beta2) = 0 + +dt_required(alg::SimpleIDSolve) = false +isdiscretealg(alg::SimpleIDSolve) = true diff --git a/lib/SimpleImplicitDiscreteSolve/src/cache.jl b/lib/SimpleImplicitDiscreteSolve/src/cache.jl new file mode 100644 index 000000000..8ecf92231 --- /dev/null +++ b/lib/SimpleImplicitDiscreteSolve/src/cache.jl @@ -0,0 +1,39 @@ +mutable struct ImplicitDiscreteState{uType, pType, tType} + u::Vector{uType} + p::Union{Nothing, Vector{pType}} + t_next::tType +end + +mutable struct SimpleIDSolveCache{uType} <: OrdinaryDiffEqMutableCache + u::uType + uprev::uType + state::ImplicitDiscreteState + prob::Union{Nothing, AbstractNonlinearProblem} +end + +function alg_cache(alg::SimpleIDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + + state = ImplicitDiscreteState(similar(u), similar(p), t) + SimpleIDSolveCache(u, uprev, state, nothing) +end + +isdiscretecache(cache::SimpleIDSolveCache) = true + +struct SimpleIDSolveConstantCache <: OrdinaryDiffEqConstantCache + prob::Union{Nothing, AbstractNonlinearProblem} +end + +function alg_cache(alg::SimpleIDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + + state = ImplicitDiscreteState(similar(u), similar(p), t) + SimpleIDSolveCache(u, uprev, state, nothing) +end + +isfsal(alg::SimpleIDSolve) = false +get_fsalfirstlast(cache::SimpleIDSolveCache, rate_prototype) = (nothing, nothing) diff --git a/lib/SimpleImplicitDiscreteSolve/src/solve.jl b/lib/SimpleImplicitDiscreteSolve/src/solve.jl new file mode 100644 index 000000000..88392faeb --- /dev/null +++ b/lib/SimpleImplicitDiscreteSolve/src/solve.jl @@ -0,0 +1,32 @@ +# Remake the nonlinear problem, then update +function perform_step!(integrator, cache::SimpleIDSolveCache, repeat_step = false) + @unpack alg, u, uprev, dt, t, f, p = integrator + @unpack state, prob = cache + state.u .= uprev + state.t_next = t + prob = remake(prob, p = state) + + u = solve(prob, SimpleNewtonRaphson()) + any(isnan, u) && (integrator.sol.retcode = SciMLBase.ReturnCode.Failure) + integrator.u = u +end + +function initialize!(integrator, cache::SimpleIDSolveCache) + cache.state.u .= integrator.u + cache.state.p .= integrator.p + cache.state.t_next = integrator.t + f = integrator.f + + _f = if isinplace(f) + (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t_next) + else + (u_next, p) -> f(u_next, p.u, p.p, p.t_next) + end + + prob = if isinplace(f) + NonlinearProblem{true}(_f, cache.state.u, cache.state) + else + NonlinearProblem{false}(_f, cache.state.u, cache.state) + end + cache.prob = prob +end From 05cd4db05885c659ad3063b223ace20fce335a52 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 26 Feb 2025 10:16:22 -0500 Subject: [PATCH 06/32] refactor --- .../src/ImplicitDiscreteSolve.jl | 32 ----------------- lib/ImplicitDiscreteSolve/src/alg_utils.jl | 19 ---------- lib/ImplicitDiscreteSolve/src/cache.jl | 35 ------------------- lib/ImplicitDiscreteSolve/src/solve.jl | 33 ----------------- .../Project.toml | 0 .../test/runtests.jl | 5 ++- 6 files changed, 2 insertions(+), 122 deletions(-) delete mode 100644 lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl delete mode 100644 lib/ImplicitDiscreteSolve/src/alg_utils.jl delete mode 100644 lib/ImplicitDiscreteSolve/src/cache.jl delete mode 100644 lib/ImplicitDiscreteSolve/src/solve.jl rename lib/{ImplicitDiscreteSolve => SimpleImplicitDiscreteSolve}/Project.toml (100%) rename lib/{ImplicitDiscreteSolve => SimpleImplicitDiscreteSolve}/test/runtests.jl (86%) diff --git a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl deleted file mode 100644 index 8fa8550d8..000000000 --- a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl +++ /dev/null @@ -1,32 +0,0 @@ -module ImplicitDiscreteSolve - -using SciMLBase: AbstractNonlinearProblem -using SciMLBase -using NonlinearSolveBase -using SymbolicIndexingInterface -using LinearAlgebra -using ADTypes -using UnPack -import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal, initialize!, perform_step!, isdiscretecache, isdiscretealg, alg_order, beta2_default, beta1_default, dt_required -import CommonSolve -import DifferentiationInterface as DI - -using Reexport -@reexport using DiffEqBase - -""" - IDSolve(alg; autodiff = true, kwargs...) - -Solver for `ImplicitDiscreteSystems`. `alg` is the NonlinearSolve algorithm that is used to solve for the next timestep at each step. -""" -struct IDSolve{algType} <: OrdinaryDiffEqAlgorithm - nlsolve::algType -end - -include("cache.jl") -include("solve.jl") -include("alg_utils.jl") - -export IDSolve - -end diff --git a/lib/ImplicitDiscreteSolve/src/alg_utils.jl b/lib/ImplicitDiscreteSolve/src/alg_utils.jl deleted file mode 100644 index 5c25c2c21..000000000 --- a/lib/ImplicitDiscreteSolve/src/alg_utils.jl +++ /dev/null @@ -1,19 +0,0 @@ -function SciMLBase.isautodifferentiable(alg::IDSolve) - true -end -function SciMLBase.allows_arbitrary_number_types(alg::IDSolve) - true -end -function SciMLBase.allowscomplex(alg::IDSolve) - true -end - -SciMLBase.isdiscrete(alg::IDSolve) = true - -isfsal(alg::IDSolve) = false -alg_order(alg::IDSolve) = 0 -beta2_default(alg::IDSolve) = 0 -beta1_default(alg::IDSolve, beta2) = 0 - -dt_required(alg::IDSolve) = false -isdiscretealg(alg::IDSolve) = true diff --git a/lib/ImplicitDiscreteSolve/src/cache.jl b/lib/ImplicitDiscreteSolve/src/cache.jl deleted file mode 100644 index 4eef79141..000000000 --- a/lib/ImplicitDiscreteSolve/src/cache.jl +++ /dev/null @@ -1,35 +0,0 @@ -mutable struct ImplicitDiscreteState{uType, pType, tType} - u::Vector{uType} - p::Union{Nothing, Vector{pType}} - t_next::tType -end - -mutable struct IDSolveCache{uType} <: OrdinaryDiffEqMutableCache - u::uType - uprev::uType - state::ImplicitDiscreteState - prob::Union{Nothing, AbstractNonlinearProblem} -end - -function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - - state = ImplicitDiscreteState(similar(u), similar(p), t) - IDSolveCache(u, uprev, state, nothing) -end - -isdiscretecache(cache::IDSolveCache) = true - -function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - - state = ImplicitDiscreteState(similar(u), similar(p), t) - IDSolveCache(u, uprev, state, nothing) -end - -isfsal(alg::IDSolve) = false -get_fsalfirstlast(cache::IDSolveCache, rate_prototype) = (nothing, nothing) diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl deleted file mode 100644 index cf76a745e..000000000 --- a/lib/ImplicitDiscreteSolve/src/solve.jl +++ /dev/null @@ -1,33 +0,0 @@ -# Remake the nonlinear problem, then update -function perform_step!(integrator, cache::IDSolveCache, repeat_step = false) - @unpack alg, u, uprev, dt, t, f, p = integrator - nlsolve = alg.nlsolve - @unpack state, prob = cache - state.u .= uprev - state.t_next = t - prob = remake(prob, p = state) - - u = solve(prob, nlsolve) - any(isnan, u) && (integrator.sol.retcode = SciMLBase.ReturnCode.Failure) - integrator.u = u -end - -function initialize!(integrator, cache::IDSolveCache) - cache.state.u .= integrator.u - cache.state.p .= integrator.p - cache.state.t_next = integrator.t - f = integrator.f - - _f = if isinplace(f) - (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t_next) - else - (u_next, p) -> f(u_next, p.u, p.p, p.t_next) - end - - prob = if isinplace(f) - NonlinearProblem{true}(_f, cache.state.u, cache.state) - else - NonlinearProblem{false}(_f, cache.state.u, cache.state) - end - cache.prob = prob -end diff --git a/lib/ImplicitDiscreteSolve/Project.toml b/lib/SimpleImplicitDiscreteSolve/Project.toml similarity index 100% rename from lib/ImplicitDiscreteSolve/Project.toml rename to lib/SimpleImplicitDiscreteSolve/Project.toml diff --git a/lib/ImplicitDiscreteSolve/test/runtests.jl b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl similarity index 86% rename from lib/ImplicitDiscreteSolve/test/runtests.jl rename to lib/SimpleImplicitDiscreteSolve/test/runtests.jl index e8bb2b601..a2397c269 100644 --- a/lib/ImplicitDiscreteSolve/test/runtests.jl +++ b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl @@ -1,8 +1,7 @@ #runtests -using ImplicitDiscreteSolve +using SimpleImplicitDiscreteSolve using OrdinaryDiffEqCore using OrdinaryDiffEqSDIRK -using SimpleNonlinearSolve # Test implicit Euler using ImplicitDiscreteProblem @testset "Implicit Discrete System" begin @@ -20,7 +19,7 @@ using SimpleNonlinearSolve tspan = (0., 0.5) idprob = ImplicitDiscreteProblem(f!, u0, tspan, []; dt = 0.01) - idsol = solve(idprob, IDSolve(SimpleNewtonRaphson())) + idsol = solve(idprob, SimpleIDSolve()) oprob = ODEProblem(lotkavolterra, u0, tspan) osol = solve(oprob, ImplicitEuler()) From d8bb4c14bc96273854b84e58cecddfcd1ab7bdd1 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 26 Feb 2025 10:58:51 -0500 Subject: [PATCH 07/32] more renaming --- lib/SimpleImplicitDiscreteSolve/Project.toml | 3 ++- .../src/SimpleImplicitDiscreteSolve.jl | 8 ++++---- lib/SimpleImplicitDiscreteSolve/src/cache.jl | 1 - 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/Project.toml b/lib/SimpleImplicitDiscreteSolve/Project.toml index 9e1096534..9273f76d8 100644 --- a/lib/SimpleImplicitDiscreteSolve/Project.toml +++ b/lib/SimpleImplicitDiscreteSolve/Project.toml @@ -1,4 +1,4 @@ -name = "ImplicitDiscreteSolve" +name = "SimpleImplicitDiscreteSolve" uuid = "3263718b-31ed-49cf-8a0f-35a466e8af96" authors = ["vyudu "] version = "0.1.0" @@ -13,6 +13,7 @@ NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" diff --git a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl index 29f61cf66..6dc741ac8 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl @@ -1,4 +1,4 @@ -module ImplicitDiscreteSolve +module SimpleImplicitDiscreteSolve using SciMLBase: AbstractNonlinearProblem using SciMLBase @@ -16,11 +16,11 @@ using Reexport @reexport using DiffEqBase """ - IDSolve(alg; autodiff = true, kwargs...) + SimpleIDSolve() -Solver for `ImplicitDiscreteSystems`. `alg` is the NonlinearSolve algorithm that is used to solve for the next timestep at each step. +Simple solver for `ImplicitDiscreteSystems`. Uses `SimpleNewtonRaphson` to solve for the next state at every timestep. """ -struct SimpleIDSolve{algType} <: OrdinaryDiffEqAlgorithm end +struct SimpleIDSolve <: OrdinaryDiffEqAlgorithm end include("cache.jl") include("solve.jl") diff --git a/lib/SimpleImplicitDiscreteSolve/src/cache.jl b/lib/SimpleImplicitDiscreteSolve/src/cache.jl index 8ecf92231..79b25e392 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/cache.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/cache.jl @@ -35,5 +35,4 @@ function alg_cache(alg::SimpleIDSolve, u, rate_prototype, ::Type{uEltypeNoUnits} SimpleIDSolveCache(u, uprev, state, nothing) end -isfsal(alg::SimpleIDSolve) = false get_fsalfirstlast(cache::SimpleIDSolveCache, rate_prototype) = (nothing, nothing) From 9780da73b186ff31f28bc083c9d34a45ba0175b7 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 1 Mar 2025 17:35:16 -0500 Subject: [PATCH 08/32] thin dependencies --- .../src/SimpleImplicitDiscreteSolve.jl | 7 ------- lib/SimpleImplicitDiscreteSolve/src/cache.jl | 4 ++-- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl index 6dc741ac8..b834a60a3 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl @@ -1,16 +1,9 @@ module SimpleImplicitDiscreteSolve -using SciMLBase: AbstractNonlinearProblem using SciMLBase -using NonlinearSolveBase using SimpleNonlinearSolve -using SymbolicIndexingInterface -using LinearAlgebra -using ADTypes using UnPack import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal, initialize!, perform_step!, isdiscretecache, isdiscretealg, alg_order, beta2_default, beta1_default, dt_required -import CommonSolve -import DifferentiationInterface as DI using Reexport @reexport using DiffEqBase diff --git a/lib/SimpleImplicitDiscreteSolve/src/cache.jl b/lib/SimpleImplicitDiscreteSolve/src/cache.jl index 79b25e392..c70b9bb53 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/cache.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/cache.jl @@ -8,7 +8,7 @@ mutable struct SimpleIDSolveCache{uType} <: OrdinaryDiffEqMutableCache u::uType uprev::uType state::ImplicitDiscreteState - prob::Union{Nothing, AbstractNonlinearProblem} + prob::Union{Nothing, SciMLBase.AbstractNonlinearProblem} end function alg_cache(alg::SimpleIDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -23,7 +23,7 @@ end isdiscretecache(cache::SimpleIDSolveCache) = true struct SimpleIDSolveConstantCache <: OrdinaryDiffEqConstantCache - prob::Union{Nothing, AbstractNonlinearProblem} + prob::Union{Nothing, SciMLBase.AbstractNonlinearProblem} end function alg_cache(alg::SimpleIDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, From 23205f5e1fe295caee417296f84560717c556feb Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 3 Mar 2025 13:20:06 -0500 Subject: [PATCH 09/32] Update Project.toml --- lib/SimpleImplicitDiscreteSolve/Project.toml | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/Project.toml b/lib/SimpleImplicitDiscreteSolve/Project.toml index 9273f76d8..a0870327e 100644 --- a/lib/SimpleImplicitDiscreteSolve/Project.toml +++ b/lib/SimpleImplicitDiscreteSolve/Project.toml @@ -4,32 +4,20 @@ authors = ["vyudu "] version = "0.1.0" [deps] -ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" -CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" -DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" -SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] -ADTypes = "1.13.0" -CommonSolve = "0.2.4" DiffEqBase = "6.164.1" -DifferentiationInterface = "0.6.42" -LinearAlgebra = "1.11.0" -NonlinearSolveBase = "1.5.0" OrdinaryDiffEqCore = "1.18.1" OrdinaryDiffEqSDIRK = "1.2.0" Reexport = "1.2.2" SciMLBase = "2.74.1" SimpleNonlinearSolve = "2.1.0" -SymbolicIndexingInterface = "0.3.37" UnPack = "1.0.2" [extras] From 918d83d769b698fa3b262392a15d524ab0d6e4fb Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 11 Mar 2025 22:48:24 -0400 Subject: [PATCH 10/32] add dep to NonlinearSolve --- Project.toml | 1 + .../test/runtests.jl | 24 ++++++++++++++++++- src/NonlinearSolve.jl | 1 + 3 files changed, 25 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f26780676..f0f781c19 100644 --- a/Project.toml +++ b/Project.toml @@ -25,6 +25,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SimpleImplicitDiscreteSolve = "3263718b-31ed-49cf-8a0f-35a466e8af96" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" diff --git a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl index a2397c269..0c8613aeb 100644 --- a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl +++ b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl @@ -4,7 +4,7 @@ using OrdinaryDiffEqCore using OrdinaryDiffEqSDIRK # Test implicit Euler using ImplicitDiscreteProblem -@testset "Implicit Discrete System" begin +@testset "Implicit Euler" begin function lotkavolterra(u, p, t) [1.5*u[1] - u[1]*u[2], -3.0*u[2] + u[1]*u[2]] end @@ -25,4 +25,26 @@ using OrdinaryDiffEqSDIRK osol = solve(oprob, ImplicitEuler()) @test isapprox(idsol[end], osol[end], atol = 0.01) + + ### free-fall + # y, dy + function ff(u, p, t) + [u[2], -9.8] + end + + function g!(resid, u_next, u, p, t) + f = ff(u_next, p, t) + resid[1] = u_next[1] - u[1] - 0.01*f[1] + resid[2] = u_next[2] - u[2] - 0.01*f[2] + nothing + end + u0 = [100., 3.] + + idprob = ImplicitDiscreteProblem(g!, u0, tspan, []; dt = 0.01) + idsol = solve(idprob, SimpleIDSolve()) + + oprob = ODEProblem(ff, u0, tspan) + osol = solve(oprob, ImplicitEuler()) + + @test isapprox(idsol[end], osol[end], atol = 0.01) end diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index c5c031b7d..e33b2669c 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -42,6 +42,7 @@ using NonlinearSolveFirstOrder: NonlinearSolveFirstOrder, GeneralizedFirstOrderA using NonlinearSolveQuasiNewton: NonlinearSolveQuasiNewton, QuasiNewtonAlgorithm using NonlinearSolveSpectralMethods: NonlinearSolveSpectralMethods, GeneralizedDFSane using SimpleNonlinearSolve: SimpleNonlinearSolve +using SimpleImplicitDiscreteSolve: SimpleIDSolve const SII = SymbolicIndexingInterface From 8ec870bbc18c92fc8691102b96857686cd06a4c8 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 12 Mar 2025 00:21:31 -0400 Subject: [PATCH 11/32] refactor: do not use similar(p) --- lib/SimpleImplicitDiscreteSolve/Project.toml | 2 ++ .../src/SimpleImplicitDiscreteSolve.jl | 1 + lib/SimpleImplicitDiscreteSolve/src/cache.jl | 6 +++--- lib/SimpleImplicitDiscreteSolve/src/solve.jl | 1 + src/NonlinearSolve.jl | 3 ++- 5 files changed, 9 insertions(+), 4 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/Project.toml b/lib/SimpleImplicitDiscreteSolve/Project.toml index a0870327e..080ea71a0 100644 --- a/lib/SimpleImplicitDiscreteSolve/Project.toml +++ b/lib/SimpleImplicitDiscreteSolve/Project.toml @@ -9,6 +9,7 @@ OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" +SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] @@ -18,6 +19,7 @@ OrdinaryDiffEqSDIRK = "1.2.0" Reexport = "1.2.2" SciMLBase = "2.74.1" SimpleNonlinearSolve = "2.1.0" +SymbolicIndexingInterface = "0.3.38" UnPack = "1.0.2" [extras] diff --git a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl index b834a60a3..09b0e0713 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl @@ -3,6 +3,7 @@ module SimpleImplicitDiscreteSolve using SciMLBase using SimpleNonlinearSolve using UnPack +using SymbolicIndexingInterface: parameter_values import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal, initialize!, perform_step!, isdiscretecache, isdiscretealg, alg_order, beta2_default, beta1_default, dt_required using Reexport diff --git a/lib/SimpleImplicitDiscreteSolve/src/cache.jl b/lib/SimpleImplicitDiscreteSolve/src/cache.jl index c70b9bb53..828ba28b9 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/cache.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/cache.jl @@ -1,6 +1,6 @@ mutable struct ImplicitDiscreteState{uType, pType, tType} u::Vector{uType} - p::Union{Nothing, Vector{pType}} + p::pType t_next::tType end @@ -16,7 +16,7 @@ function alg_cache(alg::SimpleIDSolve, u, rate_prototype, ::Type{uEltypeNoUnits} dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - state = ImplicitDiscreteState(similar(u), similar(p), t) + state = ImplicitDiscreteState(similar(u), p, t) SimpleIDSolveCache(u, uprev, state, nothing) end @@ -31,7 +31,7 @@ function alg_cache(alg::SimpleIDSolve, u, rate_prototype, ::Type{uEltypeNoUnits} dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - state = ImplicitDiscreteState(similar(u), similar(p), t) + state = ImplicitDiscreteState(similar(u), p, t) SimpleIDSolveCache(u, uprev, state, nothing) end diff --git a/lib/SimpleImplicitDiscreteSolve/src/solve.jl b/lib/SimpleImplicitDiscreteSolve/src/solve.jl index 88392faeb..8beb68957 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/solve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/solve.jl @@ -13,6 +13,7 @@ end function initialize!(integrator, cache::SimpleIDSolveCache) cache.state.u .= integrator.u + @show integrator.p cache.state.p .= integrator.p cache.state.t_next = integrator.t f = integrator.f diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index e33b2669c..206b93e58 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -110,7 +110,8 @@ end # Rexexports @reexport using SciMLBase, NonlinearSolveBase, LineSearch, ADTypes @reexport using NonlinearSolveFirstOrder, NonlinearSolveSpectralMethods, - NonlinearSolveQuasiNewton, SimpleNonlinearSolve, BracketingNonlinearSolve + NonlinearSolveQuasiNewton, SimpleNonlinearSolve, BracketingNonlinearSolve, + SimpleImplicitDiscreteSolve @reexport using LinearSolve # Poly Algorithms From d9c0782a624e36f0b4cf9822a1647a360da98b96 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 12 Mar 2025 00:38:19 -0400 Subject: [PATCH 12/32] fix: can't broadcast over MTKParameter --- .../src/SimpleImplicitDiscreteSolve.jl | 2 +- lib/SimpleImplicitDiscreteSolve/src/solve.jl | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl index 09b0e0713..f5f4f1f80 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl @@ -3,7 +3,7 @@ module SimpleImplicitDiscreteSolve using SciMLBase using SimpleNonlinearSolve using UnPack -using SymbolicIndexingInterface: parameter_values +using SymbolicIndexingInterface: parameter_symbols import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal, initialize!, perform_step!, isdiscretecache, isdiscretealg, alg_order, beta2_default, beta1_default, dt_required using Reexport diff --git a/lib/SimpleImplicitDiscreteSolve/src/solve.jl b/lib/SimpleImplicitDiscreteSolve/src/solve.jl index 8beb68957..0b37a68f6 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/solve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/solve.jl @@ -13,8 +13,7 @@ end function initialize!(integrator, cache::SimpleIDSolveCache) cache.state.u .= integrator.u - @show integrator.p - cache.state.p .= integrator.p + cache.state.p = integrator.p cache.state.t_next = integrator.t f = integrator.f From 0e56572676ca4d10bba2e419138e7910abef8049 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 12 Mar 2025 01:17:32 -0400 Subject: [PATCH 13/32] feat: add --- lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl | 15 +++++++++++++++ lib/SimpleImplicitDiscreteSolve/test/runtests.jl | 3 +++ 2 files changed, 18 insertions(+) diff --git a/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl b/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl index b4eff1bc1..c67e93ebf 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl @@ -17,3 +17,18 @@ beta1_default(alg::SimpleIDSolve, beta2) = 0 dt_required(alg::SimpleIDSolve) = false isdiscretealg(alg::SimpleIDSolve) = true + +function DiffEqBase._initialize_dae!(integrator, prob::ImplicitDiscreteProblem, + alg::DefaultInit, x::Union{Val{true}, Val{false}}) + atol = one(eltype(prob.u0)) * 1e-12 + if SciMLBase.has_initializeprob(prob.f) + _initialize_dae!(integrator, prob, + OverrideInit(atol), x) + elseif !applicable(_initialize_dae!, integrator, prob, + BrownFullBasicInit(atol), x) + error("`OrdinaryDiffEqNonlinearSolve` is not loaded, which is required for the default initialization algorithm (`BrownFullBasicInit` or `ShampineCollocationInit`). To solve this problem, either do `using OrdinaryDiffEqNonlinearSolve` or pass `initializealg = CheckInit()` to the `solve` function. This second option requires consistent `u0`.") + else + _initialize_dae!(integrator, prob, + BrownFullBasicInit(atol), x) + end +end diff --git a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl index 0c8613aeb..14dff9517 100644 --- a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl +++ b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl @@ -48,3 +48,6 @@ using OrdinaryDiffEqSDIRK @test isapprox(idsol[end], osol[end], atol = 0.01) end + +@testset "Solve respects initialization" begin +end From 5ad133c9b6916df53aa76a341a6c50703b9d104e Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 12 Mar 2025 04:01:50 -0400 Subject: [PATCH 14/32] feat: add for ImplicitDiscreteProblem --- .../src/SimpleImplicitDiscreteSolve.jl | 2 +- lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl index f5f4f1f80..151dd0993 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl @@ -4,7 +4,7 @@ using SciMLBase using SimpleNonlinearSolve using UnPack using SymbolicIndexingInterface: parameter_symbols -import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal, initialize!, perform_step!, isdiscretecache, isdiscretealg, alg_order, beta2_default, beta1_default, dt_required +import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal, initialize!, perform_step!, isdiscretecache, isdiscretealg, alg_order, beta2_default, beta1_default, dt_required, _initialize_dae!, DefaultInit, BrownFullBasicInit, OverrideInit using Reexport @reexport using DiffEqBase diff --git a/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl b/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl index c67e93ebf..53ba7666e 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl @@ -18,7 +18,7 @@ beta1_default(alg::SimpleIDSolve, beta2) = 0 dt_required(alg::SimpleIDSolve) = false isdiscretealg(alg::SimpleIDSolve) = true -function DiffEqBase._initialize_dae!(integrator, prob::ImplicitDiscreteProblem, +function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem, alg::DefaultInit, x::Union{Val{true}, Val{false}}) atol = one(eltype(prob.u0)) * 1e-12 if SciMLBase.has_initializeprob(prob.f) From 94dcaf1ef86cd36e16648d3c937425241f5466b3 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 12 Mar 2025 21:23:54 -0400 Subject: [PATCH 15/32] Add initialization --- lib/SimpleImplicitDiscreteSolve/Project.toml | 5 +- .../src/alg_utils.jl | 15 ------ lib/SimpleImplicitDiscreteSolve/src/solve.jl | 49 +++++++++++++++++++ .../test/runtests.jl | 26 ++++++++-- 4 files changed, 74 insertions(+), 21 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/Project.toml b/lib/SimpleImplicitDiscreteSolve/Project.toml index 080ea71a0..6d7b50f8c 100644 --- a/lib/SimpleImplicitDiscreteSolve/Project.toml +++ b/lib/SimpleImplicitDiscreteSolve/Project.toml @@ -20,11 +20,12 @@ Reexport = "1.2.2" SciMLBase = "2.74.1" SimpleNonlinearSolve = "2.1.0" SymbolicIndexingInterface = "0.3.38" +Test = "1.11.0" UnPack = "1.0.2" [extras] OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" -SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["OrdinaryDiffEqSDIRK", "SimpleNonlinearSolve"] +test = ["OrdinaryDiffEqSDIRK", "Test"] diff --git a/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl b/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl index 53ba7666e..b4eff1bc1 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl @@ -17,18 +17,3 @@ beta1_default(alg::SimpleIDSolve, beta2) = 0 dt_required(alg::SimpleIDSolve) = false isdiscretealg(alg::SimpleIDSolve) = true - -function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem, - alg::DefaultInit, x::Union{Val{true}, Val{false}}) - atol = one(eltype(prob.u0)) * 1e-12 - if SciMLBase.has_initializeprob(prob.f) - _initialize_dae!(integrator, prob, - OverrideInit(atol), x) - elseif !applicable(_initialize_dae!, integrator, prob, - BrownFullBasicInit(atol), x) - error("`OrdinaryDiffEqNonlinearSolve` is not loaded, which is required for the default initialization algorithm (`BrownFullBasicInit` or `ShampineCollocationInit`). To solve this problem, either do `using OrdinaryDiffEqNonlinearSolve` or pass `initializealg = CheckInit()` to the `solve` function. This second option requires consistent `u0`.") - else - _initialize_dae!(integrator, prob, - BrownFullBasicInit(atol), x) - end -end diff --git a/lib/SimpleImplicitDiscreteSolve/src/solve.jl b/lib/SimpleImplicitDiscreteSolve/src/solve.jl index 0b37a68f6..591d9ebed 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/solve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/solve.jl @@ -30,3 +30,52 @@ function initialize!(integrator, cache::SimpleIDSolveCache) end cache.prob = prob end + +function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem, + alg::DefaultInit, x::Union{Val{true}, Val{false}}) + atol = one(eltype(prob.u0)) * 1e-12 + if SciMLBase.has_initializeprob(prob.f) + _initialize_dae!(integrator, prob, + OverrideInit(atol), x) + else + @unpack u, p, t, f = integrator + initstate = ImplicitDiscreteState(u, p, t) + + _f = if isinplace(f) + (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t_next) + else + (u_next, p) -> f(u_next, p.u, p.p, p.t_next) + end + prob = NonlinearProblem{isinplace(f)}(_f, u, initstate) + sol = solve(prob, SimpleNewtonRaphson()) + integrator.u = sol + end +end + +#### TODO: Implement real algorithm +function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem, alg::BrownFullBasicInit, isinplace::Val{true}) + @unpack p, t, f = integrator + u0 = integrator.u + + nlequation! = (out, u, p) -> begin + f(out, u, u0, p, t) + end + + nlfunc = NonlinearFunction(nlequation!; jac_prototype = f.jac_prototype) + nlprob = NonlinearProblem(nlfunc, ifelse.(differential_vars, du, u), p) + nlsol = solve(nlprob, nlsolve; abstol = alg.abstol, reltol = integrator.opts.reltol) + + @. du = ifelse(differential_vars, nlsol.u, du) + @. u = ifelse(differential_vars, u, nlsol.u) + + recursivecopy!(integrator.uprev, integrator.u) + if alg_extrapolates(integrator.alg) + recursivecopy!(integrator.uprev2, integrator.uprev) + end + + if nlsol.retcode != ReturnCode.Success + integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, + ReturnCode.InitialFailure) + end + return +end diff --git a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl index 14dff9517..eea3f8289 100644 --- a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl +++ b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl @@ -1,4 +1,5 @@ #runtests +using Test using SimpleImplicitDiscreteSolve using OrdinaryDiffEqCore using OrdinaryDiffEqSDIRK @@ -24,7 +25,7 @@ using OrdinaryDiffEqSDIRK oprob = ODEProblem(lotkavolterra, u0, tspan) osol = solve(oprob, ImplicitEuler()) - @test isapprox(idsol[end], osol[end], atol = 0.01) + @test isapprox(idsol[end], osol[end], atol = 0.1) ### free-fall # y, dy @@ -38,7 +39,8 @@ using OrdinaryDiffEqSDIRK resid[2] = u_next[2] - u[2] - 0.01*f[2] nothing end - u0 = [100., 3.] + u0 = [10., 0.] + tspan = (0, 0.2) idprob = ImplicitDiscreteProblem(g!, u0, tspan, []; dt = 0.01) idsol = solve(idprob, SimpleIDSolve()) @@ -46,8 +48,24 @@ using OrdinaryDiffEqSDIRK oprob = ODEProblem(ff, u0, tspan) osol = solve(oprob, ImplicitEuler()) - @test isapprox(idsol[end], osol[end], atol = 0.01) + @test isapprox(idsol[end], osol[end], atol = 0.1) end -@testset "Solve respects initialization" begin +@testset "Solver initializes" begin + function periodic!(resid, u_next, u, p, t) + resid[1] = u_next[1] - u[1] - sin(t*π/4) + resid[2] = 16 - u_next[2]^2 - u_next[1]^2 + end + + tsteps = 15 + u0 = [1., 3.] + idprob = ImplicitDiscreteProblem(periodic!, u0, (0, tsteps), []) + integ = init(idprob, SimpleIDSolve()) + @test integ.u[1]^2 + integ.u[2]^2 ≈ 16 + + for ts in 1:tsteps + step!(integ) + @show integ.u + @test integ.u[1]^2 + integ.u[2]^2 ≈ 16 + end end From 85c33fbcbfacf7fcde628d1858b3ac16a890b4c9 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 12 Mar 2025 21:32:27 -0400 Subject: [PATCH 16/32] comment out initialize method --- lib/SimpleImplicitDiscreteSolve/src/solve.jl | 52 ++++++++++---------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/src/solve.jl b/lib/SimpleImplicitDiscreteSolve/src/solve.jl index 591d9ebed..357d511d1 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/solve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/solve.jl @@ -53,29 +53,29 @@ function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem, end #### TODO: Implement real algorithm -function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem, alg::BrownFullBasicInit, isinplace::Val{true}) - @unpack p, t, f = integrator - u0 = integrator.u - - nlequation! = (out, u, p) -> begin - f(out, u, u0, p, t) - end - - nlfunc = NonlinearFunction(nlequation!; jac_prototype = f.jac_prototype) - nlprob = NonlinearProblem(nlfunc, ifelse.(differential_vars, du, u), p) - nlsol = solve(nlprob, nlsolve; abstol = alg.abstol, reltol = integrator.opts.reltol) - - @. du = ifelse(differential_vars, nlsol.u, du) - @. u = ifelse(differential_vars, u, nlsol.u) - - recursivecopy!(integrator.uprev, integrator.u) - if alg_extrapolates(integrator.alg) - recursivecopy!(integrator.uprev2, integrator.uprev) - end - - if nlsol.retcode != ReturnCode.Success - integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, - ReturnCode.InitialFailure) - end - return -end +# function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem, alg::BrownFullBasicInit, isinplace::Val{true}) +# @unpack p, t, f = integrator +# u0 = integrator.u +# +# nlequation! = (out, u, p) -> begin +# f(out, u, u0, p, t) +# end +# +# nlfunc = NonlinearFunction(nlequation!; jac_prototype = f.jac_prototype) +# nlprob = NonlinearProblem(nlfunc, ifelse.(differential_vars, du, u), p) +# nlsol = solve(nlprob, nlsolve; abstol = alg.abstol, reltol = integrator.opts.reltol) +# +# @. du = ifelse(differential_vars, nlsol.u, du) +# @. u = ifelse(differential_vars, u, nlsol.u) +# +# recursivecopy!(integrator.uprev, integrator.u) +# if alg_extrapolates(integrator.alg) +# recursivecopy!(integrator.uprev2, integrator.uprev) +# end +# +# if nlsol.retcode != ReturnCode.Success +# integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, +# ReturnCode.InitialFailure) +# end +# return +# end From 4c54e0742169ac373b9d757c56552acd6f8544e5 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 13 Mar 2025 06:21:36 -0100 Subject: [PATCH 17/32] Update Project.toml --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index f0f781c19..f26780676 100644 --- a/Project.toml +++ b/Project.toml @@ -25,7 +25,6 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -SimpleImplicitDiscreteSolve = "3263718b-31ed-49cf-8a0f-35a466e8af96" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" From 20be3767e36d6e53755d68908db483c6224230f2 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 13 Mar 2025 06:35:06 -0100 Subject: [PATCH 18/32] Update lib/SimpleImplicitDiscreteSolve/src/solve.jl --- lib/SimpleImplicitDiscreteSolve/src/solve.jl | 27 -------------------- 1 file changed, 27 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/src/solve.jl b/lib/SimpleImplicitDiscreteSolve/src/solve.jl index 357d511d1..f6312d2cd 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/solve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/solve.jl @@ -52,30 +52,3 @@ function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem, end end -#### TODO: Implement real algorithm -# function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem, alg::BrownFullBasicInit, isinplace::Val{true}) -# @unpack p, t, f = integrator -# u0 = integrator.u -# -# nlequation! = (out, u, p) -> begin -# f(out, u, u0, p, t) -# end -# -# nlfunc = NonlinearFunction(nlequation!; jac_prototype = f.jac_prototype) -# nlprob = NonlinearProblem(nlfunc, ifelse.(differential_vars, du, u), p) -# nlsol = solve(nlprob, nlsolve; abstol = alg.abstol, reltol = integrator.opts.reltol) -# -# @. du = ifelse(differential_vars, nlsol.u, du) -# @. u = ifelse(differential_vars, u, nlsol.u) -# -# recursivecopy!(integrator.uprev, integrator.u) -# if alg_extrapolates(integrator.alg) -# recursivecopy!(integrator.uprev2, integrator.uprev) -# end -# -# if nlsol.retcode != ReturnCode.Success -# integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, -# ReturnCode.InitialFailure) -# end -# return -# end From b90176d16bb05b4ab0209afe696ad81429976cc3 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 13 Mar 2025 09:35:33 -0400 Subject: [PATCH 19/32] use u retcode --- lib/SimpleImplicitDiscreteSolve/src/solve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleImplicitDiscreteSolve/src/solve.jl b/lib/SimpleImplicitDiscreteSolve/src/solve.jl index f6312d2cd..a09754b49 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/solve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/solve.jl @@ -7,7 +7,7 @@ function perform_step!(integrator, cache::SimpleIDSolveCache, repeat_step = fals prob = remake(prob, p = state) u = solve(prob, SimpleNewtonRaphson()) - any(isnan, u) && (integrator.sol.retcode = SciMLBase.ReturnCode.Failure) + integrator.sol.retcode = u.retcode integrator.u = u end From 7faedf3d997f7fbf42547c443ef1168de3da19ff Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 24 Mar 2025 19:22:59 -0400 Subject: [PATCH 20/32] refactor: SimpleImplicitDiscreteSolve --- lib/SimpleImplicitDiscreteSolve/Project.toml | 6 -- .../src/SimpleImplicitDiscreteSolve.jl | 86 +++++++++++++++++-- .../src/alg_utils.jl | 19 ---- lib/SimpleImplicitDiscreteSolve/src/cache.jl | 38 -------- lib/SimpleImplicitDiscreteSolve/src/solve.jl | 54 ------------ .../test/runtests.jl | 16 ++-- 6 files changed, 85 insertions(+), 134 deletions(-) delete mode 100644 lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl delete mode 100644 lib/SimpleImplicitDiscreteSolve/src/cache.jl delete mode 100644 lib/SimpleImplicitDiscreteSolve/src/solve.jl diff --git a/lib/SimpleImplicitDiscreteSolve/Project.toml b/lib/SimpleImplicitDiscreteSolve/Project.toml index 6d7b50f8c..3ca1dc30f 100644 --- a/lib/SimpleImplicitDiscreteSolve/Project.toml +++ b/lib/SimpleImplicitDiscreteSolve/Project.toml @@ -5,23 +5,17 @@ version = "0.1.0" [deps] DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" -OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" -SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" -UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] DiffEqBase = "6.164.1" -OrdinaryDiffEqCore = "1.18.1" OrdinaryDiffEqSDIRK = "1.2.0" Reexport = "1.2.2" SciMLBase = "2.74.1" SimpleNonlinearSolve = "2.1.0" -SymbolicIndexingInterface = "0.3.38" Test = "1.11.0" -UnPack = "1.0.2" [extras] OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" diff --git a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl index 151dd0993..048135b40 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl @@ -2,10 +2,6 @@ module SimpleImplicitDiscreteSolve using SciMLBase using SimpleNonlinearSolve -using UnPack -using SymbolicIndexingInterface: parameter_symbols -import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal, initialize!, perform_step!, isdiscretecache, isdiscretealg, alg_order, beta2_default, beta1_default, dt_required, _initialize_dae!, DefaultInit, BrownFullBasicInit, OverrideInit - using Reexport @reexport using DiffEqBase @@ -14,11 +10,85 @@ using Reexport Simple solver for `ImplicitDiscreteSystems`. Uses `SimpleNewtonRaphson` to solve for the next state at every timestep. """ -struct SimpleIDSolve <: OrdinaryDiffEqAlgorithm end +struct SimpleIDSolve <: SciMLBase.AbstractODEAlgorithm end + +function DiffEqBase.__init(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; dt = 1) + u0 = prob.u0 + p = prob.p + f = prob.f + t = prob.tspan[1] + + nlf = isinplace(f) ? (out, u, p) -> f(out, u, u0, p, t) : (u, p) -> f(u, u0, p, t) + prob = NonlinearProblem{isinplace(f)}(nlf, u0, p) + sol = solve(prob, SimpleNewtonRaphson()) + sol, (sol.retcode != ReturnCode.Success) +end + +function DiffEqBase.solve(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; + dt = 1, + save_everystep = true, + save_start = true, + adaptive = false, + dense = false, + save_end = true, + kwargs...) + + @assert !adaptive + @assert !dense + (initsol, initfail) = DiffEqBase.__init(prob, alg; dt) + if initfail + sol = DiffEqBase.build_solution(prob, alg, prob.tspan[1], u0, k = nothing, stats = nothing, calculate_error = false) + return SciMLBase.solution_new_retcode(sol, ReturnCode.InitialFailure) + end -include("cache.jl") -include("solve.jl") -include("alg_utils.jl") + u0 = initsol.u + tspan = prob.tspan + f = prob.f + p = prob.p + t = tspan[1] + tf = prob.tspan[2] + ts = tspan[1]:dt:tspan[2] + + if save_everystep && save_start + us = Vector{typeof(u0)}(undef, length(ts)) + us[1] = u0 + elseif save_everystep + us = Vector{typeof(u0)}(undef, length(ts) - 1) + elseif save_start + us = Vector{typeof(u0)}(undef, 2) + us[1] = u0 + else + us = Vector{typeof(u0)}(undef, 1) # for interface compatibility + end + + u = u0 + convfail = false + for i in 2:length(ts) + uprev = u + t = ts[i] + nlf = isinplace(f) ? (out, u, p) -> f(out, u, uprev, p, t) : (u, p) -> f(u, uprev, p, t) + nlprob = NonlinearProblem{isinplace(f)}(nlf, uprev, p) + nlsol = solve(nlprob, SimpleNewtonRaphson()) + u = nlsol.u + save_everystep && (us[i] = u) + convfail = (nlsol.retcode != ReturnCode.Success) + + if convfail + sol = DiffEqBase.build_solution(prob, alg, ts[1:i], us[1:i], k = nothing, stats = nothing, calculate_error = false) + sol = SciMLBase.solution_new_retcode(sol, ReturnCode.ConvergenceFailure) + return sol + end + end + + !save_everystep && save_end && (us[end] = u) + sol = DiffEqBase.build_solution(prob, alg, ts, us, + k = nothing, stats = nothing, + calculate_error = false) + + DiffEqBase.has_analytic(prob.f) && + DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true, dense_errors = false) + sol +end export SimpleIDSolve diff --git a/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl b/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl deleted file mode 100644 index b4eff1bc1..000000000 --- a/lib/SimpleImplicitDiscreteSolve/src/alg_utils.jl +++ /dev/null @@ -1,19 +0,0 @@ -function SciMLBase.isautodifferentiable(alg::SimpleIDSolve) - true -end -function SciMLBase.allows_arbitrary_number_types(alg::SimpleIDSolve) - true -end -function SciMLBase.allowscomplex(alg::SimpleIDSolve) - true -end - -SciMLBase.isdiscrete(alg::SimpleIDSolve) = true - -isfsal(alg::SimpleIDSolve) = false -alg_order(alg::SimpleIDSolve) = 0 -beta2_default(alg::SimpleIDSolve) = 0 -beta1_default(alg::SimpleIDSolve, beta2) = 0 - -dt_required(alg::SimpleIDSolve) = false -isdiscretealg(alg::SimpleIDSolve) = true diff --git a/lib/SimpleImplicitDiscreteSolve/src/cache.jl b/lib/SimpleImplicitDiscreteSolve/src/cache.jl deleted file mode 100644 index 828ba28b9..000000000 --- a/lib/SimpleImplicitDiscreteSolve/src/cache.jl +++ /dev/null @@ -1,38 +0,0 @@ -mutable struct ImplicitDiscreteState{uType, pType, tType} - u::Vector{uType} - p::pType - t_next::tType -end - -mutable struct SimpleIDSolveCache{uType} <: OrdinaryDiffEqMutableCache - u::uType - uprev::uType - state::ImplicitDiscreteState - prob::Union{Nothing, SciMLBase.AbstractNonlinearProblem} -end - -function alg_cache(alg::SimpleIDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - - state = ImplicitDiscreteState(similar(u), p, t) - SimpleIDSolveCache(u, uprev, state, nothing) -end - -isdiscretecache(cache::SimpleIDSolveCache) = true - -struct SimpleIDSolveConstantCache <: OrdinaryDiffEqConstantCache - prob::Union{Nothing, SciMLBase.AbstractNonlinearProblem} -end - -function alg_cache(alg::SimpleIDSolve, u, rate_prototype, ::Type{uEltypeNoUnits}, - ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, - dt, reltol, p, calck, - ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - - state = ImplicitDiscreteState(similar(u), p, t) - SimpleIDSolveCache(u, uprev, state, nothing) -end - -get_fsalfirstlast(cache::SimpleIDSolveCache, rate_prototype) = (nothing, nothing) diff --git a/lib/SimpleImplicitDiscreteSolve/src/solve.jl b/lib/SimpleImplicitDiscreteSolve/src/solve.jl deleted file mode 100644 index a09754b49..000000000 --- a/lib/SimpleImplicitDiscreteSolve/src/solve.jl +++ /dev/null @@ -1,54 +0,0 @@ -# Remake the nonlinear problem, then update -function perform_step!(integrator, cache::SimpleIDSolveCache, repeat_step = false) - @unpack alg, u, uprev, dt, t, f, p = integrator - @unpack state, prob = cache - state.u .= uprev - state.t_next = t - prob = remake(prob, p = state) - - u = solve(prob, SimpleNewtonRaphson()) - integrator.sol.retcode = u.retcode - integrator.u = u -end - -function initialize!(integrator, cache::SimpleIDSolveCache) - cache.state.u .= integrator.u - cache.state.p = integrator.p - cache.state.t_next = integrator.t - f = integrator.f - - _f = if isinplace(f) - (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t_next) - else - (u_next, p) -> f(u_next, p.u, p.p, p.t_next) - end - - prob = if isinplace(f) - NonlinearProblem{true}(_f, cache.state.u, cache.state) - else - NonlinearProblem{false}(_f, cache.state.u, cache.state) - end - cache.prob = prob -end - -function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem, - alg::DefaultInit, x::Union{Val{true}, Val{false}}) - atol = one(eltype(prob.u0)) * 1e-12 - if SciMLBase.has_initializeprob(prob.f) - _initialize_dae!(integrator, prob, - OverrideInit(atol), x) - else - @unpack u, p, t, f = integrator - initstate = ImplicitDiscreteState(u, p, t) - - _f = if isinplace(f) - (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t_next) - else - (u_next, p) -> f(u_next, p.u, p.p, p.t_next) - end - prob = NonlinearProblem{isinplace(f)}(_f, u, initstate) - sol = solve(prob, SimpleNewtonRaphson()) - integrator.u = sol - end -end - diff --git a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl index eea3f8289..6d6ba48b3 100644 --- a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl +++ b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl @@ -1,7 +1,6 @@ #runtests using Test using SimpleImplicitDiscreteSolve -using OrdinaryDiffEqCore using OrdinaryDiffEqSDIRK # Test implicit Euler using ImplicitDiscreteProblem @@ -19,13 +18,13 @@ using OrdinaryDiffEqSDIRK u0 = [1., 1.] tspan = (0., 0.5) - idprob = ImplicitDiscreteProblem(f!, u0, tspan, []; dt = 0.01) - idsol = solve(idprob, SimpleIDSolve()) + idprob = ImplicitDiscreteProblem(f!, u0, tspan, []) + idsol = solve(idprob, SimpleIDSolve(), dt = 0.01) oprob = ODEProblem(lotkavolterra, u0, tspan) osol = solve(oprob, ImplicitEuler()) - @test isapprox(idsol[end], osol[end], atol = 0.1) + @test isapprox(idsol[end-1], osol[end], atol = 0.1) ### free-fall # y, dy @@ -40,15 +39,15 @@ using OrdinaryDiffEqSDIRK nothing end u0 = [10., 0.] - tspan = (0, 0.2) + tspan = (0, 0.5) - idprob = ImplicitDiscreteProblem(g!, u0, tspan, []; dt = 0.01) - idsol = solve(idprob, SimpleIDSolve()) + idprob = ImplicitDiscreteProblem(g!, u0, tspan, []) + idsol = solve(idprob, SimpleIDSolve(); dt = 0.01) oprob = ODEProblem(ff, u0, tspan) osol = solve(oprob, ImplicitEuler()) - @test isapprox(idsol[end], osol[end], atol = 0.1) + @test isapprox(idsol[end-1], osol[end], atol = 0.1) end @testset "Solver initializes" begin @@ -65,7 +64,6 @@ end for ts in 1:tsteps step!(integ) - @show integ.u @test integ.u[1]^2 + integ.u[2]^2 ≈ 16 end end From 8f10a85f71b78be52b734248f77168fc6eea5844 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 24 Mar 2025 19:25:25 -0400 Subject: [PATCH 21/32] fix test --- lib/SimpleImplicitDiscreteSolve/test/runtests.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl index 6d6ba48b3..d4f300595 100644 --- a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl +++ b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl @@ -59,11 +59,13 @@ end tsteps = 15 u0 = [1., 3.] idprob = ImplicitDiscreteProblem(periodic!, u0, (0, tsteps), []) - integ = init(idprob, SimpleIDSolve()) - @test integ.u[1]^2 + integ.u[2]^2 ≈ 16 + initsol, initfail = DiffEqBase.__init(idprob, SimpleIDSolve()) + @test initsol.u[1]^2 + initsol.u[2]^2 ≈ 16 + + idsol = solve(idprob, SimpleIDSolve()) for ts in 1:tsteps - step!(integ) - @test integ.u[1]^2 + integ.u[2]^2 ≈ 16 + step = idsol.u[ts] + @test step[1]^2 + step[2]^2 ≈ 16 end end From 878dc5f4643aa04ed52781b81e601b047c42cc16 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 24 Mar 2025 19:26:11 -0400 Subject: [PATCH 22/32] add Project --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index f26780676..f0f781c19 100644 --- a/Project.toml +++ b/Project.toml @@ -25,6 +25,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SimpleImplicitDiscreteSolve = "3263718b-31ed-49cf-8a0f-35a466e8af96" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" From 4f52f5ac8cafc767904134151626e5ee8e3ad670 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 21 Apr 2025 11:18:13 -0400 Subject: [PATCH 23/32] add License and CI --- .../CI_SimpleImplicitDiscreteSolve.yml | 120 ++++++++++++++++++ lib/SimpleImplicitDiscreteSolve/LICENSE | 21 +++ 2 files changed, 141 insertions(+) create mode 100644 .github/workflows/CI_SimpleImplicitDiscreteSolve.yml create mode 100644 lib/SimpleImplicitDiscreteSolve/LICENSE diff --git a/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml b/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml new file mode 100644 index 000000000..8ba5853ad --- /dev/null +++ b/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml @@ -0,0 +1,120 @@ +name: CI (SimpleNonlinearSolve) + +on: + pull_request: + branches: + - master + paths: + - "lib/SimpleImplicitDiscreteSolve/**" + - ".github/workflows/CI_SimpleImplicitDiscreteSolve.yml" + - "lib/SimpleNonlinearSolve/**" + push: + branches: + - master + +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} + +jobs: + test: + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - "1.10" + - "1" + os: + - ubuntu-latest + - macos-latest + - windows-latest + group: + - core + - adjoint + - alloc_check + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + - uses: actions/cache@v4 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - name: "Install Dependencies and Run Tests" + run: | + import Pkg + Pkg.Registry.update() + # Install packages present in subdirectories + dev_pks = Pkg.PackageSpec[] + for path in ("lib/SimpleImplicitDiscreteSolve", "lib/SimpleNonlinearSolve") + push!(dev_pks, Pkg.PackageSpec(; path)) + end + Pkg.develop(dev_pks) + Pkg.instantiate() + Pkg.test(; coverage="user") + shell: julia --color=yes --code-coverage=user --depwarn=yes --project=lib/SimpleImplicitDiscreteSolve {0} + env: + GROUP: ${{ matrix.group }} + - uses: julia-actions/julia-processcoverage@v1 + with: + directories: lib/SimpleNonlinearSolve/src, lib/SimpleNonlinearSolve/ext, lib/SimpleImplicitDiscreteSolve/src + - uses: codecov/codecov-action@v5 + with: + file: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} + verbose: true + fail_ci_if_error: false + + downgrade: + runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + version: + - "1.10" + group: + - core + - adjoint + - alloc_check + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + - uses: julia-actions/julia-downgrade-compat@v1 + with: + skip: SimpleImplicitDiscreteSolve + - name: "Install Dependencies and Run Tests" + run: | + import Pkg + Pkg.Registry.update() + # Install packages present in subdirectories + dev_pks = Pkg.PackageSpec[] + for path in ("lib/SimpleImplicitDiscreteSolve", "lib/SimpleNonlinearSolve") + push!(dev_pks, Pkg.PackageSpec(; path)) + end + Pkg.develop(dev_pks) + Pkg.instantiate() + Pkg.test(; coverage="user") + shell: julia --color=yes --code-coverage=user --depwarn=yes --project=lib/SimpleImplicitDiscreteSolve {0} + env: + GROUP: ${{ matrix.group }} + - uses: julia-actions/julia-processcoverage@v1 + with: + directories: lib/SimpleImplicitDiscreteSolve/src, lib/SimpleNonlinearSolve/ext, lib/SimpleNonlinearSolve/src + - uses: codecov/codecov-action@v5 + with: + file: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} + verbose: true + fail_ci_if_error: false diff --git a/lib/SimpleImplicitDiscreteSolve/LICENSE b/lib/SimpleImplicitDiscreteSolve/LICENSE new file mode 100644 index 000000000..abb594d1e --- /dev/null +++ b/lib/SimpleImplicitDiscreteSolve/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2024 SciML + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. From 9bff4c80b0cae14e4d72e87e5fdb7841be3faeb2 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 21 Apr 2025 11:34:05 -0400 Subject: [PATCH 24/32] fix CI script --- .github/workflows/CI_SimpleImplicitDiscreteSolve.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml b/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml index 8ba5853ad..482a0692f 100644 --- a/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml +++ b/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml @@ -1,4 +1,4 @@ -name: CI (SimpleNonlinearSolve) +name: CI (SimpleImplicitDiscreteSolve) on: pull_request: From 8663f2678bff399ff076347a9c0997f33ae5f176 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 21 Apr 2025 11:47:33 -0400 Subject: [PATCH 25/32] fix: don't load SIDS in script --- .github/workflows/CI_SimpleImplicitDiscreteSolve.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml b/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml index 482a0692f..3c89fb09e 100644 --- a/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml +++ b/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml @@ -56,7 +56,7 @@ jobs: Pkg.Registry.update() # Install packages present in subdirectories dev_pks = Pkg.PackageSpec[] - for path in ("lib/SimpleImplicitDiscreteSolve", "lib/SimpleNonlinearSolve") + for path in ("lib/SimpleNonlinearSolve",) push!(dev_pks, Pkg.PackageSpec(; path)) end Pkg.develop(dev_pks) @@ -100,7 +100,7 @@ jobs: Pkg.Registry.update() # Install packages present in subdirectories dev_pks = Pkg.PackageSpec[] - for path in ("lib/SimpleImplicitDiscreteSolve", "lib/SimpleNonlinearSolve") + for path in ("lib/SimpleNonlinearSolve",) push!(dev_pks, Pkg.PackageSpec(; path)) end Pkg.develop(dev_pks) From fe465cc7dc7550d59517d3dbd1e4bfef10844f04 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 21 Apr 2025 11:52:08 -0400 Subject: [PATCH 26/32] fix test compat --- lib/SimpleImplicitDiscreteSolve/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleImplicitDiscreteSolve/Project.toml b/lib/SimpleImplicitDiscreteSolve/Project.toml index 3ca1dc30f..0574f11c2 100644 --- a/lib/SimpleImplicitDiscreteSolve/Project.toml +++ b/lib/SimpleImplicitDiscreteSolve/Project.toml @@ -15,7 +15,7 @@ OrdinaryDiffEqSDIRK = "1.2.0" Reexport = "1.2.2" SciMLBase = "2.74.1" SimpleNonlinearSolve = "2.1.0" -Test = "1.11.0" +Test = "1.10" [extras] OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" From f6b25793beed2f5d1a61d29c1db9b6da77a56fff Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 21 Apr 2025 12:07:34 -0400 Subject: [PATCH 27/32] change SimpleImplicitDiscreteSolve UUID --- .github/workflows/CI_SimpleImplicitDiscreteSolve.yml | 4 ++-- Project.toml | 2 +- lib/SimpleImplicitDiscreteSolve/Project.toml | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml b/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml index 3c89fb09e..4f05ce758 100644 --- a/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml +++ b/.github/workflows/CI_SimpleImplicitDiscreteSolve.yml @@ -67,7 +67,7 @@ jobs: GROUP: ${{ matrix.group }} - uses: julia-actions/julia-processcoverage@v1 with: - directories: lib/SimpleNonlinearSolve/src, lib/SimpleNonlinearSolve/ext, lib/SimpleImplicitDiscreteSolve/src + directories: lib/SimpleNonlinearSolve/src,lib/SimpleNonlinearSolve/ext,lib/SimpleImplicitDiscreteSolve/src - uses: codecov/codecov-action@v5 with: file: lcov.info @@ -111,7 +111,7 @@ jobs: GROUP: ${{ matrix.group }} - uses: julia-actions/julia-processcoverage@v1 with: - directories: lib/SimpleImplicitDiscreteSolve/src, lib/SimpleNonlinearSolve/ext, lib/SimpleNonlinearSolve/src + directories: lib/SimpleImplicitDiscreteSolve/src,lib/SimpleNonlinearSolve/ext,lib/SimpleNonlinearSolve/src - uses: codecov/codecov-action@v5 with: file: lcov.info diff --git a/Project.toml b/Project.toml index f0f781c19..2c59c9cd5 100644 --- a/Project.toml +++ b/Project.toml @@ -25,7 +25,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -SimpleImplicitDiscreteSolve = "3263718b-31ed-49cf-8a0f-35a466e8af96" +SimpleImplicitDiscreteSolve = "8b67ef88-54bd-43ff-aca0-8be02588656a" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" diff --git a/lib/SimpleImplicitDiscreteSolve/Project.toml b/lib/SimpleImplicitDiscreteSolve/Project.toml index 0574f11c2..7f21e0e05 100644 --- a/lib/SimpleImplicitDiscreteSolve/Project.toml +++ b/lib/SimpleImplicitDiscreteSolve/Project.toml @@ -1,5 +1,5 @@ name = "SimpleImplicitDiscreteSolve" -uuid = "3263718b-31ed-49cf-8a0f-35a466e8af96" +uuid = "8b67ef88-54bd-43ff-aca0-8be02588656a" authors = ["vyudu "] version = "0.1.0" From aa126f77390357a2e00d02ef9cade73c7540fa60 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 21 Apr 2025 12:16:30 -0400 Subject: [PATCH 28/32] format --- .../src/SimpleImplicitDiscreteSolve.jl | 29 +++++++++------- .../test/runtests.jl | 34 +++++++++---------- 2 files changed, 33 insertions(+), 30 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl index 048135b40..0a417b3fe 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl @@ -24,20 +24,20 @@ function DiffEqBase.__init(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; dt sol, (sol.retcode != ReturnCode.Success) end -function DiffEqBase.solve(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; - dt = 1, - save_everystep = true, - save_start = true, - adaptive = false, - dense = false, - save_end = true, - kwargs...) - +function DiffEqBase.solve(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; + dt = 1, + save_everystep = true, + save_start = true, + adaptive = false, + dense = false, + save_end = true, + kwargs...) @assert !adaptive @assert !dense (initsol, initfail) = DiffEqBase.__init(prob, alg; dt) if initfail - sol = DiffEqBase.build_solution(prob, alg, prob.tspan[1], u0, k = nothing, stats = nothing, calculate_error = false) + sol = DiffEqBase.build_solution(prob, alg, prob.tspan[1], u0, k = nothing, + stats = nothing, calculate_error = false) return SciMLBase.solution_new_retcode(sol, ReturnCode.InitialFailure) end @@ -66,7 +66,8 @@ function DiffEqBase.solve(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; for i in 2:length(ts) uprev = u t = ts[i] - nlf = isinplace(f) ? (out, u, p) -> f(out, u, uprev, p, t) : (u, p) -> f(u, uprev, p, t) + nlf = isinplace(f) ? (out, u, p) -> f(out, u, uprev, p, t) : + (u, p) -> f(u, uprev, p, t) nlprob = NonlinearProblem{isinplace(f)}(nlf, uprev, p) nlsol = solve(nlprob, SimpleNewtonRaphson()) u = nlsol.u @@ -74,7 +75,8 @@ function DiffEqBase.solve(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; convfail = (nlsol.retcode != ReturnCode.Success) if convfail - sol = DiffEqBase.build_solution(prob, alg, ts[1:i], us[1:i], k = nothing, stats = nothing, calculate_error = false) + sol = DiffEqBase.build_solution(prob, alg, ts[1:i], us[1:i], k = nothing, + stats = nothing, calculate_error = false) sol = SciMLBase.solution_new_retcode(sol, ReturnCode.ConvergenceFailure) return sol end @@ -86,7 +88,8 @@ function DiffEqBase.solve(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; calculate_error = false) DiffEqBase.has_analytic(prob.f) && - DiffEqBase.calculate_solution_errors!(sol; timeseries_errors = true, dense_errors = false) + DiffEqBase.calculate_solution_errors!( + sol; timeseries_errors = true, dense_errors = false) sol end diff --git a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl index d4f300595..0607fbe69 100644 --- a/lib/SimpleImplicitDiscreteSolve/test/runtests.jl +++ b/lib/SimpleImplicitDiscreteSolve/test/runtests.jl @@ -6,17 +6,17 @@ using OrdinaryDiffEqSDIRK # Test implicit Euler using ImplicitDiscreteProblem @testset "Implicit Euler" begin function lotkavolterra(u, p, t) - [1.5*u[1] - u[1]*u[2], -3.0*u[2] + u[1]*u[2]] + [1.5 * u[1] - u[1] * u[2], -3.0 * u[2] + u[1] * u[2]] end function f!(resid, u_next, u, p, t) lv = lotkavolterra(u_next, p, t) - resid[1] = u_next[1] - u[1] - 0.01*lv[1] - resid[2] = u_next[2] - u[2] - 0.01*lv[2] + resid[1] = u_next[1] - u[1] - 0.01 * lv[1] + resid[2] = u_next[2] - u[2] - 0.01 * lv[2] nothing end - u0 = [1., 1.] - tspan = (0., 0.5) + u0 = [1.0, 1.0] + tspan = (0.0, 0.5) idprob = ImplicitDiscreteProblem(f!, u0, tspan, []) idsol = solve(idprob, SimpleIDSolve(), dt = 0.01) @@ -24,21 +24,21 @@ using OrdinaryDiffEqSDIRK oprob = ODEProblem(lotkavolterra, u0, tspan) osol = solve(oprob, ImplicitEuler()) - @test isapprox(idsol[end-1], osol[end], atol = 0.1) + @test isapprox(idsol[end - 1], osol[end], atol = 0.1) ### free-fall # y, dy - function ff(u, p, t) + function ff(u, p, t) [u[2], -9.8] end - function g!(resid, u_next, u, p, t) - f = ff(u_next, p, t) - resid[1] = u_next[1] - u[1] - 0.01*f[1] - resid[2] = u_next[2] - u[2] - 0.01*f[2] + function g!(resid, u_next, u, p, t) + f = ff(u_next, p, t) + resid[1] = u_next[1] - u[1] - 0.01 * f[1] + resid[2] = u_next[2] - u[2] - 0.01 * f[2] nothing end - u0 = [10., 0.] + u0 = [10.0, 0.0] tspan = (0, 0.5) idprob = ImplicitDiscreteProblem(g!, u0, tspan, []) @@ -47,17 +47,17 @@ using OrdinaryDiffEqSDIRK oprob = ODEProblem(ff, u0, tspan) osol = solve(oprob, ImplicitEuler()) - @test isapprox(idsol[end-1], osol[end], atol = 0.1) + @test isapprox(idsol[end - 1], osol[end], atol = 0.1) end @testset "Solver initializes" begin - function periodic!(resid, u_next, u, p, t) - resid[1] = u_next[1] - u[1] - sin(t*π/4) - resid[2] = 16 - u_next[2]^2 - u_next[1]^2 + function periodic!(resid, u_next, u, p, t) + resid[1] = u_next[1] - u[1] - sin(t * π / 4) + resid[2] = 16 - u_next[2]^2 - u_next[1]^2 end tsteps = 15 - u0 = [1., 3.] + u0 = [1.0, 3.0] idprob = ImplicitDiscreteProblem(periodic!, u0, (0, tsteps), []) initsol, initfail = DiffEqBase.__init(idprob, SimpleIDSolve()) @test initsol.u[1]^2 + initsol.u[2]^2 ≈ 16 From 6725fc878472bef304e6dff5e3925035e468cd19 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 21 Apr 2025 12:22:15 -0400 Subject: [PATCH 29/32] Update src/NonlinearSolve.jl --- src/NonlinearSolve.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 206b93e58..405b64714 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -42,7 +42,6 @@ using NonlinearSolveFirstOrder: NonlinearSolveFirstOrder, GeneralizedFirstOrderA using NonlinearSolveQuasiNewton: NonlinearSolveQuasiNewton, QuasiNewtonAlgorithm using NonlinearSolveSpectralMethods: NonlinearSolveSpectralMethods, GeneralizedDFSane using SimpleNonlinearSolve: SimpleNonlinearSolve -using SimpleImplicitDiscreteSolve: SimpleIDSolve const SII = SymbolicIndexingInterface From d532041bc4c3f7e1685b74741531dd54ca56e0eb Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 21 Apr 2025 12:22:33 -0400 Subject: [PATCH 30/32] Update src/NonlinearSolve.jl --- src/NonlinearSolve.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 405b64714..c5c031b7d 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -109,8 +109,7 @@ end # Rexexports @reexport using SciMLBase, NonlinearSolveBase, LineSearch, ADTypes @reexport using NonlinearSolveFirstOrder, NonlinearSolveSpectralMethods, - NonlinearSolveQuasiNewton, SimpleNonlinearSolve, BracketingNonlinearSolve, - SimpleImplicitDiscreteSolve + NonlinearSolveQuasiNewton, SimpleNonlinearSolve, BracketingNonlinearSolve @reexport using LinearSolve # Poly Algorithms From 6cba7b722693f8a7ec22be28aba360bd35a55e7f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 21 Apr 2025 12:24:01 -0400 Subject: [PATCH 31/32] Update Project.toml --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 2c59c9cd5..f26780676 100644 --- a/Project.toml +++ b/Project.toml @@ -25,7 +25,6 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -SimpleImplicitDiscreteSolve = "8b67ef88-54bd-43ff-aca0-8be02588656a" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" From 5d119bcaef3502988c14ab38f8926fd6579efb1e Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 21 Apr 2025 13:06:43 -0400 Subject: [PATCH 32/32] add StaticArrays --- lib/SimpleImplicitDiscreteSolve/Project.toml | 2 ++ .../src/SimpleImplicitDiscreteSolve.jl | 16 +++++++--------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/lib/SimpleImplicitDiscreteSolve/Project.toml b/lib/SimpleImplicitDiscreteSolve/Project.toml index 7f21e0e05..a065b7f27 100644 --- a/lib/SimpleImplicitDiscreteSolve/Project.toml +++ b/lib/SimpleImplicitDiscreteSolve/Project.toml @@ -8,6 +8,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] DiffEqBase = "6.164.1" @@ -15,6 +16,7 @@ OrdinaryDiffEqSDIRK = "1.2.0" Reexport = "1.2.2" SciMLBase = "2.74.1" SimpleNonlinearSolve = "2.1.0" +StaticArrays = "1.9.13" Test = "1.10" [extras] diff --git a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl index 0a417b3fe..8e3e8a523 100644 --- a/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl +++ b/lib/SimpleImplicitDiscreteSolve/src/SimpleImplicitDiscreteSolve.jl @@ -3,6 +3,7 @@ module SimpleImplicitDiscreteSolve using SciMLBase using SimpleNonlinearSolve using Reexport +using StaticArrays @reexport using DiffEqBase """ @@ -49,16 +50,13 @@ function DiffEqBase.solve(prob::ImplicitDiscreteProblem, alg::SimpleIDSolve; tf = prob.tspan[2] ts = tspan[1]:dt:tspan[2] - if save_everystep && save_start - us = Vector{typeof(u0)}(undef, length(ts)) - us[1] = u0 - elseif save_everystep - us = Vector{typeof(u0)}(undef, length(ts) - 1) - elseif save_start - us = Vector{typeof(u0)}(undef, 2) + l = save_everystep ? length(ts) - 1 : 1 + save_start && (l = l + 1) + u0type = typeof(u0) + us = u0type <: StaticArray ? MVector{l, u0type}(undef) : Vector{u0type}(undef, l) + + if save_start us[1] = u0 - else - us = Vector{typeof(u0)}(undef, 1) # for interface compatibility end u = u0