From add7059c8a52df7401aee9af3f7e6d8833b59db4 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 28 Oct 2024 17:14:30 -0400 Subject: [PATCH 01/24] make Observable definition comply with new record_from_solution definition --- ext/MTKBifurcationKitExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl index a0bd3bfc02..e5b57efc5a 100644 --- a/ext/MTKBifurcationKitExt.jl +++ b/ext/MTKBifurcationKitExt.jl @@ -59,7 +59,7 @@ struct ObservableRecordFromSolution{S, T} end end # Functor function that computes the value. -function (orfs::ObservableRecordFromSolution)(x, p) +function (orfs::ObservableRecordFromSolution)(x, p; k...) # Updates the state values (in subs_vals). for state_idx in 1:(orfs.state_end_idxs) orfs.subs_vals[state_idx] = orfs.subs_vals[state_idx][1] => x[state_idx] From eb681c10403031e3a1ee28c96f9d4b2defc3e389 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 28 Oct 2024 17:26:25 -0400 Subject: [PATCH 02/24] add test --- test/extensions/bifurcationkit.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/extensions/bifurcationkit.jl b/test/extensions/bifurcationkit.jl index b8c95dca99..629edf46a6 100644 --- a/test/extensions/bifurcationkit.jl +++ b/test/extensions/bifurcationkit.jl @@ -162,4 +162,12 @@ let bf = bifurcationdiagram(bp, PALC(), 2, opts_br) @test bf.γ.specialpoint[1].param≈0.1 atol=1e-4 rtol=1e-4 + + # Test with plot variable as observable + pvar = ModelingToolkit.get_var_to_name(fol)[:RHS] + bp = BifurcationProblem(fol, u0, par, bif_par; plot_var = pvar) + opts_br = ContinuationPar(p_min = -1.0, + p_max = 1.0) + bf = bifurcationdiagram(bp, PALC(), 2, opts_br) + @test bf.γ.specialpoint[1].param≈0.1 atol=1e-4 rtol=1e-4 end From 9733460668a662b41fd4cf87102d08521c1e4f92 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 22 Nov 2024 13:25:30 -0500 Subject: [PATCH 03/24] init --- src/systems/diffeqs/abstractodesystem.jl | 110 +++++++++++++++++++++++ 1 file changed, 110 insertions(+) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 8d9e0b5381..2579ebafcf 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -466,6 +466,116 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, initializeprobpmap = initializeprobpmap) end +""" +```julia +SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + jac = true, sparse = true, + simplify = false, + kwargs...) where {iip} +``` + +Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and +`ps` are used to set the order of the dependent variable and parameter vectors, +respectively. +""" +function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) + BVProblem{true}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem(sys::AbstractODESystem, + u0map::StaticArray, + args...; + kwargs...) + BVProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) +end + +function SciMLBase.BVProblem{true}(sys::AbstractODESystem, args...; kwargs...) + BVProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) + BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], + tspan = get_tspan(sys), + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + jac = true, sparse = true, + sparsity = true, + callback = nothing, + check_length = true, + warn_initialize_determined = true, + eval_expression = false, + eval_module = @__MODULE__, + kwargs...) where {iip, specialize} + + if !iscomplete(sys) + error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") + end + + f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; + t = tspan !== nothing ? tspan[1] : tspan, + check_length, warn_initialize_determined, eval_expression, eval_module, jac, sparse, sparsity, kwargs...) + + # if jac + # jac_gen = generate_jacobian(sys, dvs, ps; + # simplify = simplify, sparse = sparse, + # expression = Val{true}, + # expression_module = eval_module, + # checkbounds = checkbounds, kwargs...) + # jac_oop, jac_iip = eval_or_rgf.(jac_gen; eval_expression, eval_module) + + # _jac(u, p, t) = jac_oop(u, p, t) + # _jac(J, u, p, t) = jac_iip(J, u, p, t) + # _jac(u, p::Tuple{Vararg{Number}}, t) = jac_oop(u, p, t) + # _jac(J, u, p::Tuple{Vararg{Number}}, t) = jac_iip(J, u, p, t) + # _jac(u, p::Tuple, t) = jac_oop(u, p..., t) + # _jac(J, u, p::Tuple, t) = jac_iip(J, u, p..., t) + # _jac(u, p::MTKParameters, t) = jac_oop(u, p..., t) + # _jac(J, u, p::MTKParameters, t) = jac_iip(J, u, p..., t) + # else + # _jac = nothing + # end + + # jac_prototype = if sparse + # uElType = u0 === nothing ? Float64 : eltype(u0) + # if jac + # similar(calculate_jacobian(sys, sparse = sparse), uElType) + # else + # similar(jacobian_sparsity(sys), uElType) + # end + # else + # nothing + # end + + # f.jac = _jac + # f.jac_prototype = jac_prototype + # f.sparsity = jac ? jacobian_sparsity(sys) : nothing + + cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) + + kwargs = filter_kwargs(kwargs) + + kwargs1 = (;) + if cbs !== nothing + kwargs1 = merge(kwargs1, (callback = cbs,)) + end + + # Define the boundary conditions + bc = if iip + (residual, u, p, t) -> (residual = u[1] - u0) + else + (u, p, t) -> (u[1] - u0) + end + + return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) +end + +get_callback(prob::BVProblem) = prob.kwargs[:callback] + """ ```julia DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys), From b3da8137a9288b9c3eeff209a541dd41a5f97524 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sun, 1 Dec 2024 16:18:18 -0500 Subject: [PATCH 04/24] up --- src/systems/diffeqs/abstractodesystem.jl | 64 +++++++++--------------- 1 file changed, 23 insertions(+), 41 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 7a8bf56562..bb4cbbb6c3 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -469,6 +469,17 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, initializeprobpmap = initializeprobpmap) end +""" +```julia +SciMLBase.BVPFunction{iip}(sys::AbstractODESystem, u0map, tspan, + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + jac = true, sparse = true, + simplify = false, + kwargs...) where {iip} +``` +""" + """ ```julia SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, @@ -481,7 +492,7 @@ SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and `ps` are used to set the order of the dependent variable and parameter vectors, -respectively. +respectively. `u0` should be either the initial condition, a vector of values `u(t_i)` for collocation methods, or a function returning one or the other. """ function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) BVProblem{true}(sys, args...; kwargs...) @@ -502,12 +513,13 @@ function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) end +# figure out what's going on when we try to set `sparse`? + function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], tspan = get_tspan(sys), parammap = DiffEqBase.NullParameters(); version = nothing, tgrad = false, jac = true, sparse = true, - sparsity = true, callback = nothing, check_length = true, warn_initialize_determined = true, @@ -521,57 +533,27 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; t = tspan !== nothing ? tspan[1] : tspan, - check_length, warn_initialize_determined, eval_expression, eval_module, jac, sparse, sparsity, kwargs...) - - # if jac - # jac_gen = generate_jacobian(sys, dvs, ps; - # simplify = simplify, sparse = sparse, - # expression = Val{true}, - # expression_module = eval_module, - # checkbounds = checkbounds, kwargs...) - # jac_oop, jac_iip = eval_or_rgf.(jac_gen; eval_expression, eval_module) - - # _jac(u, p, t) = jac_oop(u, p, t) - # _jac(J, u, p, t) = jac_iip(J, u, p, t) - # _jac(u, p::Tuple{Vararg{Number}}, t) = jac_oop(u, p, t) - # _jac(J, u, p::Tuple{Vararg{Number}}, t) = jac_iip(J, u, p, t) - # _jac(u, p::Tuple, t) = jac_oop(u, p..., t) - # _jac(J, u, p::Tuple, t) = jac_iip(J, u, p..., t) - # _jac(u, p::MTKParameters, t) = jac_oop(u, p..., t) - # _jac(J, u, p::MTKParameters, t) = jac_iip(J, u, p..., t) - # else - # _jac = nothing - # end - - # jac_prototype = if sparse - # uElType = u0 === nothing ? Float64 : eltype(u0) - # if jac - # similar(calculate_jacobian(sys, sparse = sparse), uElType) - # else - # similar(jacobian_sparsity(sys), uElType) - # end - # else - # nothing - # end - - # f.jac = _jac - # f.jac_prototype = jac_prototype - # f.sparsity = jac ? jacobian_sparsity(sys) : nothing + check_length, warn_initialize_determined, eval_expression, eval_module, jac, kwargs...) cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) - kwargs = filter_kwargs(kwargs) kwargs1 = (;) if cbs !== nothing kwargs1 = merge(kwargs1, (callback = cbs,)) end + + # Construct initial conditions + _u0 = prepare_initial_state(u0) + __u0 = if _u0 isa Function + _u0(t_i) + end # Define the boundary conditions bc = if iip - (residual, u, p, t) -> (residual = u[1] - u0) + (residual, u, p, t) -> (residual = u[1] - __u0) else - (u, p, t) -> (u[1] - u0) + (u, p, t) -> (u[1] - __u0) end return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) From 4affeac4b340a06c770373b498ec6b0e94c25a06 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sun, 1 Dec 2024 17:35:05 -0500 Subject: [PATCH 05/24] up --- src/systems/diffeqs/abstractodesystem.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 263a75d153..33bddece30 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -545,9 +545,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] # Construct initial conditions _u0 = prepare_initial_state(u0) - __u0 = if _u0 isa Function - _u0(t_i) - end + __u0 = _u0 isa Function ? _u0(tspan[1]) : _u0 # Define the boundary conditions bc = if iip From a3429ea2b7d9e67898023eac1599e71c3d1b7bec Mon Sep 17 00:00:00 2001 From: vyudu Date: Sun, 1 Dec 2024 17:42:39 -0500 Subject: [PATCH 06/24] up --- src/systems/diffeqs/abstractodesystem.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 33bddece30..84f40a01f7 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -513,8 +513,6 @@ function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) end -# figure out what's going on when we try to set `sparse`? - function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], tspan = get_tspan(sys), parammap = DiffEqBase.NullParameters(); @@ -544,14 +542,13 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] end # Construct initial conditions - _u0 = prepare_initial_state(u0) - __u0 = _u0 isa Function ? _u0(tspan[1]) : _u0 + _u0 = u0 isa Function ? u0(tspan[1]) : u0 # Define the boundary conditions bc = if iip - (residual, u, p, t) -> (residual = u[1] - __u0) + (residual, u, p, t) -> (residual = u[1] - _u0) else - (u, p, t) -> (u[1] - __u0) + (u, p, t) -> (u[1] - _u0) end return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) From f751fbb51c49e0b35ee84a12c54ec2de919660a5 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sun, 1 Dec 2024 21:22:08 -0500 Subject: [PATCH 07/24] up --- src/systems/diffeqs/abstractodesystem.jl | 18 ++------ test/bvproblem.jl | 56 ++++++++++++++++++++++++ 2 files changed, 59 insertions(+), 15 deletions(-) create mode 100644 test/bvproblem.jl diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 84f40a01f7..112419dea8 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -469,17 +469,6 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, initializeprobpmap = initializeprobpmap) end -""" -```julia -SciMLBase.BVPFunction{iip}(sys::AbstractODESystem, u0map, tspan, - parammap = DiffEqBase.NullParameters(); - version = nothing, tgrad = false, - jac = true, sparse = true, - simplify = false, - kwargs...) where {iip} -``` -""" - """ ```julia SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, @@ -492,7 +481,7 @@ SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and `ps` are used to set the order of the dependent variable and parameter vectors, -respectively. `u0` should be either the initial condition, a vector of values `u(t_i)` for collocation methods, or a function returning one or the other. +respectively. `u0map` should be used to specify the initial condition, or be a function returning an initial condition. """ function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) BVProblem{true}(sys, args...; kwargs...) @@ -517,7 +506,6 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] tspan = get_tspan(sys), parammap = DiffEqBase.NullParameters(); version = nothing, tgrad = false, - jac = true, sparse = true, callback = nothing, check_length = true, warn_initialize_determined = true, @@ -531,7 +519,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; t = tspan !== nothing ? tspan[1] : tspan, - check_length, warn_initialize_determined, eval_expression, eval_module, jac, kwargs...) + check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) kwargs = filter_kwargs(kwargs) @@ -546,7 +534,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] # Define the boundary conditions bc = if iip - (residual, u, p, t) -> (residual = u[1] - _u0) + (residual, u, p, t) -> residual .= u[1] - _u0 else (u, p, t) -> (u[1] - _u0) end diff --git a/test/bvproblem.jl b/test/bvproblem.jl new file mode 100644 index 0000000000..90d41a96ad --- /dev/null +++ b/test/bvproblem.jl @@ -0,0 +1,56 @@ +using BoundaryValueDiffEq, OrdinaryDiffEq +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D + +@parameters σ = 10. ρ = 28 β = 8/3 +@variables x(t) = 1 y(t) = 0 z(t) = 0 + +eqs = [D(x) ~ σ*(y-x), + D(y) ~ x*(ρ-z)-y, + D(z) ~ x*y - β*z] + +u0map = [:x => 1., :y => 0., :z => 0.] +parammap = [:ρ => 28., :β => 8/3, :σ => 10.] +tspan = (0., 10.) + +@mtkbuild lorenz = ODESystem(eqs, t) + +bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lorenz, u0map, tspan, parammap) +sol = solve(bvp, MIRK4(), dt = 0.1); + +bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lorenz, u0map, tspan, parammap) +sol2 = solve(bvp, MIRK4(), dt = 0.1); + +op = ODEProblem(lorenz, u0map, tspan, parammap) +osol = solve(op) + +@test sol.u[end] ≈ osol.u[end] +@test sol2.u[end] ≈ osol.u[end] +@test sol.u[1] == [1., 0., 0.] +@test sol2.u[1] == [1., 0., 0.] + +### Testing on pendulum + +@parameters g = 9.81 L = 1. +@variables θ(t) = π/2 + +eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] + +@mtkbuild pend = ODESystem(eqs, t) + +u0map = [θ => π/2, D(θ) => π/2] +parammap = [:L => 2.] +tspan = (0., 10.) + +bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) +sol = solve(bvp, MIRK4(), dt = 0.05); + +bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) +sol2 = solve(bvp2, MIRK4(), dt = 0.05); + +osol = solve(pend) + +@test sol.u[end] ≈ osol.u[end] +@test sol.u[1] == [π/2, π/2] +@test sol2.u[end] ≈ osol.u[end] +@test sol2.u[1] == [π/2, π/2] From a9fdfd6a3a114c4df28e4781c6b667990c01bafc Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 3 Dec 2024 05:04:45 -0500 Subject: [PATCH 08/24] up --- src/systems/diffeqs/abstractodesystem.jl | 6 ++-- test/bvproblem.jl | 44 ++++++++++++------------ 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 112419dea8..b795bb981d 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -529,12 +529,12 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] kwargs1 = merge(kwargs1, (callback = cbs,)) end - # Construct initial conditions + # Construct initial conditions. _u0 = u0 isa Function ? u0(tspan[1]) : u0 - # Define the boundary conditions + # Define the boundary conditions. bc = if iip - (residual, u, p, t) -> residual .= u[1] - _u0 + (residual, u, p, t) -> (residual .= u[1] - _u0) else (u, p, t) -> (u[1] - _u0) end diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 90d41a96ad..4865c867b3 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -2,32 +2,31 @@ using BoundaryValueDiffEq, OrdinaryDiffEq using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D -@parameters σ = 10. ρ = 28 β = 8/3 -@variables x(t) = 1 y(t) = 0 z(t) = 0 +@parameters α = 7.5 β = 4. γ = 8. δ = 5. +@variables x(t) = 1. y(t) = 2. -eqs = [D(x) ~ σ*(y-x), - D(y) ~ x*(ρ-z)-y, - D(z) ~ x*y - β*z] +eqs = [D(x) ~ α*x - β*x*y, + D(y) ~ -γ*y + δ*x*y] -u0map = [:x => 1., :y => 0., :z => 0.] -parammap = [:ρ => 28., :β => 8/3, :σ => 10.] +u0map = [:x => 1., :y => 2.] +parammap = [:α => 7.5, :β => 4, :γ => 8., :δ => 5.] tspan = (0., 10.) -@mtkbuild lorenz = ODESystem(eqs, t) +@mtkbuild lotkavolterra = ODESystem(eqs, t) -bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lorenz, u0map, tspan, parammap) -sol = solve(bvp, MIRK4(), dt = 0.1); +bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap) +sol = solve(bvp, MIRK4(), dt = 0.01); -bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lorenz, u0map, tspan, parammap) -sol2 = solve(bvp, MIRK4(), dt = 0.1); +bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap) +sol2 = solve(bvp, MIRK4(), dt = 0.01); -op = ODEProblem(lorenz, u0map, tspan, parammap) -osol = solve(op) +op = ODEProblem(lotkavolterra, u0map, tspan, parammap) +osol = solve(op, Vern9()) -@test sol.u[end] ≈ osol.u[end] -@test sol2.u[end] ≈ osol.u[end] -@test sol.u[1] == [1., 0., 0.] -@test sol2.u[1] == [1., 0., 0.] +@test isapprox(sol.u[end],osol.u[end]; atol = 0.001) +@test isapprox(sol2.u[end],osol.u[end]; atol = 0.001) +@test sol.u[1] == [1., 2.] +@test sol2.u[1] == [1., 2.] ### Testing on pendulum @@ -39,16 +38,17 @@ eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] @mtkbuild pend = ODESystem(eqs, t) u0map = [θ => π/2, D(θ) => π/2] -parammap = [:L => 2.] +parammap = [:L => 2., :g => 9.81] tspan = (0., 10.) bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) -sol = solve(bvp, MIRK4(), dt = 0.05); +sol = solve(bvp, MIRK4(), dt = 0.01); bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) -sol2 = solve(bvp2, MIRK4(), dt = 0.05); +sol2 = solve(bvp2, MIRK4(), dt = 0.01); -osol = solve(pend) +op = ODEProblem(pend, u0map, tspan, parammap) +osol = solve(op, Vern9()) @test sol.u[end] ≈ osol.u[end] @test sol.u[1] == [π/2, π/2] From a9f210691c1f38ef0f575eb89a939763b81c9c1b Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 3 Dec 2024 19:11:35 -0500 Subject: [PATCH 09/24] up --- src/systems/diffeqs/abstractodesystem.jl | 4 ++-- test/bvproblem.jl | 8 ++++---- test/runtests.jl | 1 + 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index b795bb981d..8dac19f296 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -534,7 +534,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] # Define the boundary conditions. bc = if iip - (residual, u, p, t) -> (residual .= u[1] - _u0) + (residual, u, p, t) -> (residual .= u[1] .- _u0) else (u, p, t) -> (u[1] - _u0) end @@ -542,7 +542,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) end -get_callback(prob::BVProblem) = prob.kwargs[:callback] +get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") """ ```julia diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 4865c867b3..7235638cf0 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -23,8 +23,8 @@ sol2 = solve(bvp, MIRK4(), dt = 0.01); op = ODEProblem(lotkavolterra, u0map, tspan, parammap) osol = solve(op, Vern9()) -@test isapprox(sol.u[end],osol.u[end]; atol = 0.001) -@test isapprox(sol2.u[end],osol.u[end]; atol = 0.001) +@test isapprox(sol.u[end],osol.u[end]; atol = 0.01) +@test isapprox(sol2.u[end],osol.u[end]; atol = 0.01) @test sol.u[1] == [1., 2.] @test sol2.u[1] == [1., 2.] @@ -50,7 +50,7 @@ sol2 = solve(bvp2, MIRK4(), dt = 0.01); op = ODEProblem(pend, u0map, tspan, parammap) osol = solve(op, Vern9()) -@test sol.u[end] ≈ osol.u[end] +@test isapprox(sol.u[end], osol.u[end]; atol = 0.01) @test sol.u[1] == [π/2, π/2] -@test sol2.u[end] ≈ osol.u[end] +@test isapprox(sol2.u[end], osol.u[end]; atol = 0.01) @test sol2.u[1] == [π/2, π/2] diff --git a/test/runtests.jl b/test/runtests.jl index 44846eed57..eaa87e6407 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -80,6 +80,7 @@ end @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl") @safetestset "PDE Construction Test" include("pde.jl") @safetestset "JumpSystem Test" include("jumpsystem.jl") + @safetestset "BVProblem Test" include("bvproblem.jl") @safetestset "print_tree" include("print_tree.jl") @safetestset "Constraints Test" include("constraints.jl") end From 18fdd5f79702627e30ce0acc971677a5cd1d303f Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 13 Dec 2024 10:04:59 +0900 Subject: [PATCH 10/24] up --- Project.toml | 6 ++++-- test/bvproblem.jl | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 093c632c42..19063a0342 100644 --- a/Project.toml +++ b/Project.toml @@ -72,14 +72,15 @@ MTKBifurcationKitExt = "BifurcationKit" MTKChainRulesCoreExt = "ChainRulesCore" MTKDeepDiffsExt = "DeepDiffs" MTKHomotopyContinuationExt = "HomotopyContinuation" -MTKLabelledArraysExt = "LabelledArrays" MTKInfiniteOptExt = "InfiniteOpt" +MTKLabelledArraysExt = "LabelledArrays" [compat] AbstractTrees = "0.3, 0.4" ArrayInterface = "6, 7" BifurcationKit = "0.4" BlockArrays = "1.1" +BoundaryValueDiffEq = "5.12.0" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" @@ -145,6 +146,7 @@ julia = "1.9" [extras] AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" @@ -174,4 +176,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET"] +test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEq", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET"] diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 7235638cf0..7787bd9c3e 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -44,7 +44,7 @@ tspan = (0., 10.) bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) sol = solve(bvp, MIRK4(), dt = 0.01); -bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) +bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) sol2 = solve(bvp2, MIRK4(), dt = 0.01); op = ODEProblem(pend, u0map, tspan, parammap) From 9d65a3345ade530654a64ace2069ff24e8c67875 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 16 Dec 2024 23:43:35 +0800 Subject: [PATCH 11/24] fixing create_array --- Project.toml | 3 ++ src/systems/diffeqs/abstractodesystem.jl | 10 ++++- test/bvproblem.jl | 53 +++++++++++++++--------- 3 files changed, 45 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 19063a0342..ab854b80bd 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,9 @@ version = "9.54.0" AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" +BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" +BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a" +BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 8dac19f296..3a502d3183 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -539,11 +539,19 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] (u, p, t) -> (u[1] - _u0) end - return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) + return BVProblem{iip}(f, bc, _u0, tspan, p; kwargs1..., kwargs...) end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") +@inline function create_array(::Type{Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where dims + [elems...] +end + +@inline function create_array(::Type{Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where dims + T[elems...] +end + """ ```julia DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys), diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 7787bd9c3e..7cdb7e1837 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -2,6 +2,8 @@ using BoundaryValueDiffEq, OrdinaryDiffEq using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D +solvers = [MIRK4, RadauIIa5, LobattoIIIa3] + @parameters α = 7.5 β = 4. γ = 8. δ = 5. @variables x(t) = 1. y(t) = 2. @@ -13,20 +15,26 @@ parammap = [:α => 7.5, :β => 4, :γ => 8., :δ => 5.] tspan = (0., 10.) @mtkbuild lotkavolterra = ODESystem(eqs, t) +op = ODEProblem(lotkavolterra, u0map, tspan, parammap) +osol = solve(op, Vern9()) -bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap) -sol = solve(bvp, MIRK4(), dt = 0.01); +bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true) -bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap) -sol2 = solve(bvp, MIRK4(), dt = 0.01); +for solver in solvers + println("$solver") + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [1., 2.] +end -op = ODEProblem(lotkavolterra, u0map, tspan, parammap) -osol = solve(op, Vern9()) +# Test out of place +bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true) -@test isapprox(sol.u[end],osol.u[end]; atol = 0.01) -@test isapprox(sol2.u[end],osol.u[end]; atol = 0.01) -@test sol.u[1] == [1., 2.] -@test sol2.u[1] == [1., 2.] +for solver in solvers + sol = solve(bvp2, solver(), dt = 0.01) + @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) + @test sol.u[1] == [1., 2.] +end ### Testing on pendulum @@ -38,19 +46,24 @@ eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] @mtkbuild pend = ODESystem(eqs, t) u0map = [θ => π/2, D(θ) => π/2] -parammap = [:L => 2., :g => 9.81] +parammap = [:L => 1., :g => 9.81] tspan = (0., 10.) +op = ODEProblem(pend, u0map, tspan, parammap) +osol = solve(op, Vern9()) + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) -sol = solve(bvp, MIRK4(), dt = 0.01); +for solver in solvers + sol = solve(bvp2, solver(), dt = 0.01) + @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) + @test sol.u[1] == [π/2, π/2] +end +# Test out-of-place bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) -sol2 = solve(bvp2, MIRK4(), dt = 0.01); - -op = ODEProblem(pend, u0map, tspan, parammap) -osol = solve(op, Vern9()) -@test isapprox(sol.u[end], osol.u[end]; atol = 0.01) -@test sol.u[1] == [π/2, π/2] -@test isapprox(sol2.u[end], osol.u[end]; atol = 0.01) -@test sol2.u[1] == [π/2, π/2] +for solver in solvers + sol = solve(bvp2, solver(), dt = 0.01) + @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) + @test sol.u[1] == [π/2, π/2] +end From 999ec308d58e32c51941743a6ba5a8f096c56ac2 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 16 Dec 2024 23:44:24 +0800 Subject: [PATCH 12/24] revert Project.toml --- Project.toml | 3 --- 1 file changed, 3 deletions(-) diff --git a/Project.toml b/Project.toml index ab854b80bd..19063a0342 100644 --- a/Project.toml +++ b/Project.toml @@ -7,9 +7,6 @@ version = "9.54.0" AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" -BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" -BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a" -BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" From 9226ad687ff06cbc1dbbdcaa16ba3df85ef32d4e Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 16 Dec 2024 23:56:04 +0800 Subject: [PATCH 13/24] Up --- test/bvproblem.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 7cdb7e1837..86e3722eec 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -21,7 +21,6 @@ osol = solve(op, Vern9()) bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true) for solver in solvers - println("$solver") sol = solve(bvp, solver(), dt = 0.01) @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) @test sol.u[1] == [1., 2.] @@ -47,15 +46,15 @@ eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] u0map = [θ => π/2, D(θ) => π/2] parammap = [:L => 1., :g => 9.81] -tspan = (0., 10.) +tspan = (0., 6.) op = ODEProblem(pend, u0map, tspan, parammap) osol = solve(op, Vern9()) bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) for solver in solvers - sol = solve(bvp2, solver(), dt = 0.01) - @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) @test sol.u[1] == [π/2, π/2] end From 67d8164c591b74a9b985ffbcaf1ef248fb0efaaa Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 17 Dec 2024 00:09:52 +0800 Subject: [PATCH 14/24] formatting --- src/systems/diffeqs/abstractodesystem.jl | 11 ++++--- test/bvproblem.jl | 42 +++++++++++++----------- 2 files changed, 28 insertions(+), 25 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 87ee6e823d..06c83073bf 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -512,7 +512,6 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] eval_expression = false, eval_module = @__MODULE__, kwargs...) where {iip, specialize} - if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") end @@ -528,12 +527,12 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] if cbs !== nothing kwargs1 = merge(kwargs1, (callback = cbs,)) end - + # Construct initial conditions. _u0 = u0 isa Function ? u0(tspan[1]) : u0 # Define the boundary conditions. - bc = if iip + bc = if iip (residual, u, p, t) -> (residual .= u[1] .- _u0) else (u, p, t) -> (u[1] - _u0) @@ -544,11 +543,13 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") -@inline function create_array(::Type{Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where dims +@inline function create_array(::Type{Base.ReinterpretArray}, ::Nothing, + ::Val{1}, ::Val{dims}, elems...) where {dims} [elems...] end -@inline function create_array(::Type{Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where dims +@inline function create_array( + ::Type{Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} T[elems...] end diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 86e3722eec..1072874917 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -4,49 +4,51 @@ using ModelingToolkit: t_nounits as t, D_nounits as D solvers = [MIRK4, RadauIIa5, LobattoIIIa3] -@parameters α = 7.5 β = 4. γ = 8. δ = 5. -@variables x(t) = 1. y(t) = 2. +@parameters α=7.5 β=4.0 γ=8.0 δ=5.0 +@variables x(t)=1.0 y(t)=2.0 -eqs = [D(x) ~ α*x - β*x*y, - D(y) ~ -γ*y + δ*x*y] +eqs = [D(x) ~ α * x - β * x * y, + D(y) ~ -γ * y + δ * x * y] -u0map = [:x => 1., :y => 2.] -parammap = [:α => 7.5, :β => 4, :γ => 8., :δ => 5.] -tspan = (0., 10.) +u0map = [:x => 1.0, :y => 2.0] +parammap = [:α => 7.5, :β => 4, :γ => 8.0, :δ => 5.0] +tspan = (0.0, 10.0) @mtkbuild lotkavolterra = ODESystem(eqs, t) op = ODEProblem(lotkavolterra, u0map, tspan, parammap) osol = solve(op, Vern9()) -bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true) +bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}( + lotkavolterra, u0map, tspan, parammap; eval_expression = true) for solver in solvers sol = solve(bvp, solver(), dt = 0.01) @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - @test sol.u[1] == [1., 2.] + @test sol.u[1] == [1.0, 2.0] end # Test out of place -bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true) +bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}( + lotkavolterra, u0map, tspan, parammap; eval_expression = true) for solver in solvers sol = solve(bvp2, solver(), dt = 0.01) - @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) - @test sol.u[1] == [1., 2.] + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [1.0, 2.0] end ### Testing on pendulum -@parameters g = 9.81 L = 1. -@variables θ(t) = π/2 +@parameters g=9.81 L=1.0 +@variables θ(t) = π / 2 eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] @mtkbuild pend = ODESystem(eqs, t) -u0map = [θ => π/2, D(θ) => π/2] -parammap = [:L => 1., :g => 9.81] -tspan = (0., 6.) +u0map = [θ => π / 2, D(θ) => π / 2] +parammap = [:L => 1.0, :g => 9.81] +tspan = (0.0, 6.0) op = ODEProblem(pend, u0map, tspan, parammap) osol = solve(op, Vern9()) @@ -55,7 +57,7 @@ bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, pa for solver in solvers sol = solve(bvp, solver(), dt = 0.01) @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - @test sol.u[1] == [π/2, π/2] + @test sol.u[1] == [π / 2, π / 2] end # Test out-of-place @@ -63,6 +65,6 @@ bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, for solver in solvers sol = solve(bvp2, solver(), dt = 0.01) - @test isapprox(sol.u[end],osol.u[end]; atol = 0.01) - @test sol.u[1] == [π/2, π/2] + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [π / 2, π / 2] end From 25988f3bc6b6d66b008b6a40341f75e4547e574a Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 17 Dec 2024 10:52:37 +0800 Subject: [PATCH 15/24] up --- src/systems/diffeqs/abstractodesystem.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 06c83073bf..c84f5ff5be 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -543,13 +543,13 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") -@inline function create_array(::Type{Base.ReinterpretArray}, ::Nothing, +@inline function SciML.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where {dims} [elems...] end -@inline function create_array( - ::Type{Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} +@inline function SciML.Code.create_array( + ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} T[elems...] end From bb28d4fe2a0d753221104186fe42808c7657f5d6 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 17 Dec 2024 10:53:22 +0800 Subject: [PATCH 16/24] up --- src/systems/diffeqs/abstractodesystem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index c84f5ff5be..eac2df16aa 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -543,12 +543,12 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") -@inline function SciML.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, +@inline function SciMLBase.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where {dims} [elems...] end -@inline function SciML.Code.create_array( +@inline function SciMLBase.Code.create_array( ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} T[elems...] end From b2bf7c05532bdf68d93e566d07c5e7444473dd7a Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 17 Dec 2024 11:00:30 +0800 Subject: [PATCH 17/24] fix --- src/systems/diffeqs/abstractodesystem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index eac2df16aa..3d143d49fb 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -543,12 +543,12 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") -@inline function SciMLBase.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, +@inline function SymbolicUtils.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where {dims} [elems...] end -@inline function SciMLBase.Code.create_array( +@inline function SymbolicUtils.Code.create_array( ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} T[elems...] end From 3751c2a92c282593ef52b67b727c7db635ea677e Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 20 Dec 2024 23:58:28 +0900 Subject: [PATCH 18/24] up --- src/systems/diffeqs/abstractodesystem.jl | 165 +++++++++++------------ test/bvproblem.jl | 4 +- 2 files changed, 83 insertions(+), 86 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 3d143d49fb..d8ff71c324 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -469,90 +469,6 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, initializeprobpmap = initializeprobpmap) end -""" -```julia -SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, - parammap = DiffEqBase.NullParameters(); - version = nothing, tgrad = false, - jac = true, sparse = true, - simplify = false, - kwargs...) where {iip} -``` - -Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and -`ps` are used to set the order of the dependent variable and parameter vectors, -respectively. `u0map` should be used to specify the initial condition, or be a function returning an initial condition. -""" -function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) - BVProblem{true}(sys, args...; kwargs...) -end - -function SciMLBase.BVProblem(sys::AbstractODESystem, - u0map::StaticArray, - args...; - kwargs...) - BVProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) -end - -function SciMLBase.BVProblem{true}(sys::AbstractODESystem, args...; kwargs...) - BVProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) -end - -function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) - BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) -end - -function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], - tspan = get_tspan(sys), - parammap = DiffEqBase.NullParameters(); - version = nothing, tgrad = false, - callback = nothing, - check_length = true, - warn_initialize_determined = true, - eval_expression = false, - eval_module = @__MODULE__, - kwargs...) where {iip, specialize} - if !iscomplete(sys) - error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") - end - - f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; - t = tspan !== nothing ? tspan[1] : tspan, - check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) - - cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) - kwargs = filter_kwargs(kwargs) - - kwargs1 = (;) - if cbs !== nothing - kwargs1 = merge(kwargs1, (callback = cbs,)) - end - - # Construct initial conditions. - _u0 = u0 isa Function ? u0(tspan[1]) : u0 - - # Define the boundary conditions. - bc = if iip - (residual, u, p, t) -> (residual .= u[1] .- _u0) - else - (u, p, t) -> (u[1] - _u0) - end - - return BVProblem{iip}(f, bc, _u0, tspan, p; kwargs1..., kwargs...) -end - -get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") - -@inline function SymbolicUtils.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, - ::Val{1}, ::Val{dims}, elems...) where {dims} - [elems...] -end - -@inline function SymbolicUtils.Code.create_array( - ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} - T[elems...] -end - """ ```julia DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys), @@ -943,6 +859,87 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map = end get_callback(prob::ODEProblem) = prob.kwargs[:callback] +""" +```julia +SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + jac = true, sparse = true, + simplify = false, + kwargs...) where {iip} +``` + +Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and +`ps` are used to set the order of the dependent variable and parameter vectors, +respectively. `u0map` should be used to specify the initial condition. +""" +function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) + BVProblem{true}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem(sys::AbstractODESystem, + u0map::StaticArray, + args...; + kwargs...) + BVProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) +end + +function SciMLBase.BVProblem{true}(sys::AbstractODESystem, args...; kwargs...) + BVProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) + BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], + tspan = get_tspan(sys), + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + callback = nothing, + check_length = true, + warn_initialize_determined = true, + eval_expression = false, + eval_module = @__MODULE__, + kwargs...) where {iip, specialize} + if !iscomplete(sys) + error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") + end + + f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; + t = tspan !== nothing ? tspan[1] : tspan, + check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) + + cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) + kwargs = filter_kwargs(kwargs) + + kwargs1 = (;) + if cbs !== nothing + kwargs1 = merge(kwargs1, (callback = cbs,)) + end + + # Define the boundary conditions. + bc = if iip + (residual, u, p, t) -> (residual .= u[1] .- u0) + else + (u, p, t) -> (u[1] - u0) + end + + return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) +end + +get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") + +@inline function SymbolicUtils.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, + ::Val{1}, ::Val{dims}, elems...) where {dims} + [elems...] +end + +@inline function SymbolicUtils.Code.create_array( + ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} + T[elems...] +end + """ ```julia DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan, diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 1072874917..c5a302147d 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -10,8 +10,8 @@ solvers = [MIRK4, RadauIIa5, LobattoIIIa3] eqs = [D(x) ~ α * x - β * x * y, D(y) ~ -γ * y + δ * x * y] -u0map = [:x => 1.0, :y => 2.0] -parammap = [:α => 7.5, :β => 4, :γ => 8.0, :δ => 5.0] +u0map = [x => 1.0, y => 2.0] +parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0] tspan = (0.0, 10.0) @mtkbuild lotkavolterra = ODESystem(eqs, t) From ef1f089cbd493272e2b9297c2189aa2a6029d72f Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 8 Jan 2025 15:12:12 -0500 Subject: [PATCH 19/24] up --- src/systems/diffeqs/abstractodesystem.jl | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index d8ff71c324..87ab83f10b 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -930,16 +930,6 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") -@inline function SymbolicUtils.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing, - ::Val{1}, ::Val{dims}, elems...) where {dims} - [elems...] -end - -@inline function SymbolicUtils.Code.create_array( - ::Type{<:Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where {dims} - T[elems...] -end - """ ```julia DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan, From 2a25200bcb6a1ad57c4329b6383ba844eb73c2f7 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 9 Jan 2025 17:36:38 -0500 Subject: [PATCH 20/24] extend BVProblem for constraint equations --- src/systems/diffeqs/abstractodesystem.jl | 74 +++++++++++++++++-- test/bvproblem.jl | 90 ++++++++++++++++++++++++ 2 files changed, 160 insertions(+), 4 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index f5c1c288d7..74b1bf7596 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -873,7 +873,7 @@ function SciMLBase.BVProblem(sys::AbstractODESystem, kwargs...) BVProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) end - +o function SciMLBase.BVProblem{true}(sys::AbstractODESystem, args...; kwargs...) BVProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) end @@ -908,11 +908,32 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] kwargs1 = merge(kwargs1, (callback = cbs,)) end + # Handle algebraic equations + stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))]) + pidxmap = Dict([v => i for (i, v) in enumerate(parameters(sys))]) + ns = length(stmap) + ne = length(get_alg_eqs(sys)) + # Define the boundary conditions. - bc = if iip - (residual, u, p, t) -> (residual .= u[1] .- u0) + bc = if has_alg_eqs(sys) + if iip + (residual,u,p,t) -> begin + residual[1:ns] .= u[1] .- u0 + residual[ns+1:ns+ne] .= sub_u_p_into_symeq.(get_alg_eqs(sys)) + end + else + (u,p,t) -> begin + resid = vcat(u[1] - u0, sub_u_p_into_symeq.(get_alg_eqs(sys))) + end + end else - (u, p, t) -> (u[1] - u0) + if iip + (residual,u,p,t) -> begin + residual .= u[1] .- u0 + end + else + (u,p,t) -> (u[1] - u0) + end end return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) @@ -920,6 +941,51 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") +# Helper to create the dictionary that will substitute numeric values for u, p into the algebraic equations in the ODESystem. Used to construct the boundary condition function. +# Take a system with variables x,y, parameters g +# +# 1 + x + y → 1 + u[1][1] + u[1][2] +# x(0.5) → u(0.5)[1] +# x(0.5)*g(0.5) → u(0.5)[1]*p[1] + +function sub_u_p_into_symeq(eq, u, p, stidxmap, pidxmap) + iv = ModelingToolkit.get_iv(sys) + eq = Symbolics.unwrap(eq) + + stmap = Dict([st => u[1][i] for st => i in stidxmap]) + pmap = Dict([pa => p[i] for pa => i in pidxmap]) + eq = Symbolics.substitute(eq, merge(stmap, pmap)) + + csyms = [] + # Find most nested calls, substitute those first. + while !isempty(find_callable_syms!(csyms, eq)) + for sym in csyms + t = arguments(sym)[1] + x = operation(sym) + + if isparameter(x) + eq = Symbolics.substitute(eq, Dict(x(t) => p[pidxmap[x(iv)]])) + elseif isvariable(x) + eq = Symbolics.substitute(eq, Dict(x(t) => u(val)[stidxmap[x(iv)]])) + end + end + empty!(csyms) + end + eq +end + +function find_callable_syms!(csyms, ex) + ex = Symbolics.unwrap(ex) + + if iscall(ex) + operation(ex) isa Symbolic && (arguments(ex)[1] isa Symbolic) && push!(csyms, ex) # only add leaf nodes + for arg in arguments(ex) + find_callable_syms!(csyms, arg) + end + end + csyms +end + """ ```julia DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan, diff --git a/test/bvproblem.jl b/test/bvproblem.jl index c5a302147d..2d5535325a 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -2,6 +2,7 @@ using BoundaryValueDiffEq, OrdinaryDiffEq using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D +### Test Collocation solvers on simple problems solvers = [MIRK4, RadauIIa5, LobattoIIIa3] @parameters α=7.5 β=4.0 γ=8.0 δ=5.0 @@ -68,3 +69,92 @@ for solver in solvers @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) @test sol.u[1] == [π / 2, π / 2] end + +################################################### +### TESTING ODESystem with Constraint Equations ### +################################################### + +# Cartesian pendulum from the docs. Testing that initialization is satisfied. +let + @parameters g + @variables x(t) y(t) [state_priority = 10] λ(t) + eqs = [D(D(x)) ~ λ * x + D(D(y)) ~ λ * y - g + x^2 + y^2 ~ 1] + @mtkbuild pend = ODESystem(eqs, t) + + tspan = (0.0, 1.5) + u0map = [x => 1, y => 0] + parammap = [g => 1] + guesses = [λ => 1] + + prob = ODEProblem(pend, u0map, tspan, pmap; guesses) + sol = solve(prob, Rodas5P()) + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses) + + for solver in solvers + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + conditions = getfield.(equations(pend)[3:end], :rhs) + @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 + end + + bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) + for solver in solvers + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + conditions = getfield.(equations(pend)[3:end], :rhs) + @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 + end +end + +# Adding a midpoint boundary condition. +let + @parameters g + @variables x(..) y(t) [state_priority = 10] λ(t) + eqs = [D(D(x(t))) ~ λ * x(t) + D(D(y)) ~ λ * y - g + x(t)^2 + y^2 ~ 1 + x(0.5) ~ 1] + @mtkbuild pend = ODESystem(eqs, t) + + tspan = (0.0, 1.5) + u0map = [x(t) => 0.6, y => 0.8] + parammap = [g => 1] + guesses = [λ => 1] + + prob = ODEProblem(pend, u0map, tspan, pmap; guesses, check_length = false) + sol = solve(prob, Rodas5P()) + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses = guesses, check_length = false) + + for solver in solvers + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + conditions = getfield.(equations(pend)[3:end], :rhs) + @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 + @test sol.u[1] == [π / 2, π / 2] + end + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses) + + for solver in solvers + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + conditions = getfield.(equations(pend)[3:end], :rhs) + @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 + end + + bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) + for solver in solvers + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + conditions = getfield.(equations(pend)[3:end], :rhs) + @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 + end +end + +# Testing a more complicated case with multiple constraints. +let +end From 50504abbcf0149d50bd6e858ae1c5f368f8d2835 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 11 Jan 2025 13:57:59 -0500 Subject: [PATCH 21/24] adding tests --- Project.toml | 6 +- src/systems/diffeqs/abstractodesystem.jl | 140 ++++++++++--- test/bvproblem.jl | 242 ++++++++++++++--------- 3 files changed, 261 insertions(+), 127 deletions(-) diff --git a/Project.toml b/Project.toml index 98fd119f0c..b0bd4381b6 100644 --- a/Project.toml +++ b/Project.toml @@ -82,6 +82,7 @@ ArrayInterface = "6, 7" BifurcationKit = "0.4" BlockArrays = "1.1" BoundaryValueDiffEq = "5.12.0" +BoundaryValueDiffEqAscher = "1.1.0" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" @@ -140,8 +141,8 @@ SimpleNonlinearSolve = "0.1.0, 1, 2" SparseArrays = "1" SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2" StaticArrays = "0.10, 0.11, 0.12, 1.0" -StochasticDiffEq = "6.72.1" StochasticDelayDiffEq = "1.8.1" +StochasticDiffEq = "6.72.1" SymbolicIndexingInterface = "0.3.36" SymbolicUtils = "3.7" Symbolics = "6.19" @@ -154,6 +155,7 @@ julia = "1.9" AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" +BoundaryValueDiffEqAscher = "7227322d-7511-4e07-9247-ad6ff830280e" ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" @@ -185,4 +187,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve"] +test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEq", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve"] diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 74b1bf7596..3c45b1b2f8 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -853,15 +853,47 @@ get_callback(prob::ODEProblem) = prob.kwargs[:callback] ```julia SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, parammap = DiffEqBase.NullParameters(); + constraints = nothing, guesses = nothing, version = nothing, tgrad = false, jac = true, sparse = true, simplify = false, kwargs...) where {iip} ``` -Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and +Create a boundary value problem from the [`ODESystem`](@ref). The arguments `dvs` and `ps` are used to set the order of the dependent variable and parameter vectors, -respectively. `u0map` should be used to specify the initial condition. +respectively. `u0map` is used to specify fixed initial values for the states. + +Every variable must have either an initial guess supplied using `guesses` or +a fixed initial value specified using `u0map`. + +`constraints` are used to specify boundary conditions to the ODESystem in the +form of equations. These values should specify values that state variables should +take at specific points, as in `x(0.5) ~ 1`). More general constraints that +should hold over the entire solution, such as `x(t)^2 + y(t)^2`, should be +specified as one of the equations used to build the `ODESystem`. Below is an example. + +```julia + @parameters g + @variables x(..) y(t) [state_priority = 10] λ(t) + eqs = [D(D(x(t))) ~ λ * x(t) + D(D(y)) ~ λ * y - g + x(t)^2 + y^2 ~ 1] + @mtkbuild pend = ODESystem(eqs, t) + + tspan = (0.0, 1.5) + u0map = [x(t) => 0.6, y => 0.8] + parammap = [g => 1] + guesses = [λ => 1] + constraints = [x(0.5) ~ 1] + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) +``` + +If no `constraints` are specified, the problem will be treated as an initial value problem. + +If the `ODESystem` has algebraic equations like `x(t)^2 + y(t)^2`, the resulting +`BVProblem` must be solved using BVDAE solvers, such as Ascher. """ function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) BVProblem{true}(sys, args...; kwargs...) @@ -873,7 +905,7 @@ function SciMLBase.BVProblem(sys::AbstractODESystem, kwargs...) BVProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) end -o + function SciMLBase.BVProblem{true}(sys::AbstractODESystem, args...; kwargs...) BVProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) end @@ -885,6 +917,7 @@ end function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], tspan = get_tspan(sys), parammap = DiffEqBase.NullParameters(); + constraints = nothing, guesses = nothing, version = nothing, tgrad = false, callback = nothing, check_length = true, @@ -892,38 +925,63 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] eval_expression = false, eval_module = @__MODULE__, kwargs...) where {iip, specialize} + if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") end + !isnothing(callbacks) && error("BVP solvers do not support callbacks.") - f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; - t = tspan !== nothing ? tspan[1] : tspan, - check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) + iv = get_iv(sys) + constraintsts = nothing + constraintps = nothing + sts = unknowns(sys) + ps = parameters(sys) - cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...) - kwargs = filter_kwargs(kwargs) + if !isnothing(constraints) + constraints isa Equation || + constraints isa Vector{Equation} || + error("Constraints must be specified as an equation or a vector of equations.") - kwargs1 = (;) - if cbs !== nothing - kwargs1 = merge(kwargs1, (callback = cbs,)) + (length(constraints) + length(u0map) > length(sts)) && + error("The BVProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) cannot exceed the total number of states.") + + constraintsts = OrderedSet() + constraintps = OrderedSet() + + for eq in constraints + collect_vars!(constraintsts, constraintps, eq, iv) + validate_constraint_syms(eq, constraintsts, constraintps, Set(sts), Set(ps), iv) + empty!(constraintsts) + empty!(constraintps) + end end - # Handle algebraic equations - stidxmap = Dict([v => i for (i, v) in enumerate(unknowns(sys))]) - pidxmap = Dict([v => i for (i, v) in enumerate(parameters(sys))]) - ns = length(stmap) - ne = length(get_alg_eqs(sys)) + f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; + t = tspan !== nothing ? tspan[1] : tspan, guesses, + check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) + + stidxmap = Dict([v => i for (i, v) in enumerate(sts)]) + pidxmap = Dict([v => i for (i, v) in enumerate(ps)]) + + # Indices of states that have initial constraints. + u0i = has_alg_eqs(sys) ? collect(1:length(sts)) : [stidxmap[k] for k in keys(u0map)] + ni = length(u0i) - # Define the boundary conditions. - bc = if has_alg_eqs(sys) + bc = if !isnothing(constraints) + ne = length(constraints) if iip (residual,u,p,t) -> begin - residual[1:ns] .= u[1] .- u0 - residual[ns+1:ns+ne] .= sub_u_p_into_symeq.(get_alg_eqs(sys)) + residual[1:ni] .= u[1][u0i] .- u0[u0i] + residual[ni+1:ni+ne] .= map(constraints) do cons + sub_u_p_into_symeq(cons.rhs - cons.lhs, u, p, stidxmap, pidxmap, iv, tspan) + end end else (u,p,t) -> begin - resid = vcat(u[1] - u0, sub_u_p_into_symeq.(get_alg_eqs(sys))) + consresid = map(constraints) do cons + sub_u_p_into_symeq(cons.rhs-cons.lhs, u, p, stidxmap, pidxmap, iv, tspan) + end + resid = vcat(u[1][u0i] - u0[u0i], consresid) end end else @@ -941,32 +999,54 @@ end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") -# Helper to create the dictionary that will substitute numeric values for u, p into the algebraic equations in the ODESystem. Used to construct the boundary condition function. +# Validate that all the variables in the BVP constraints are well-formed states or parameters. +function validate_constraint_syms(eq, constraintsts, constraintps, sts, ps, iv) + ModelingToolkit.check_variables(constraintsts) + ModelingToolkit.check_parameters(constraintps) + + for var in constraintsts + if arguments(var) == iv + var ∈ sts || error("Constraint equation $eq contains a variable $var that is not a variable of the ODESystem.") + error("Constraint equation $eq contains a variable $var that does not have a specified argument. Such equations should be specified as algebraic equations to the ODESystem rather than a boundary constraints.") + else + operation(var)(iv) ∈ sts || error("Constraint equation $eq contains a variable $(operation(var)) that is not a variable of the ODESystem.") + end + end + + for var in constraintps + if !iscall(var) + var ∈ ps || error("Constraint equation $eq contains a parameter $var that is not a parameter of the ODESystem.") + else + operation(var) ∈ ps || error("Constraint equations contain a parameter $var that is not a parameter of the ODESystem.") + end + end +end + +# Helper to substitute numeric values for u, p into the algebraic equations in the ODESystem. Used to construct the boundary condition function. # Take a system with variables x,y, parameters g # -# 1 + x + y → 1 + u[1][1] + u[1][2] +# 1 + x(0) + y(0) → 1 + u[1][1] + u[1][2] # x(0.5) → u(0.5)[1] # x(0.5)*g(0.5) → u(0.5)[1]*p[1] - -function sub_u_p_into_symeq(eq, u, p, stidxmap, pidxmap) - iv = ModelingToolkit.get_iv(sys) +function sub_u_p_into_symeq(eq, u, p, stidxmap, pidxmap, iv, tspan) eq = Symbolics.unwrap(eq) - stmap = Dict([st => u[1][i] for st => i in stidxmap]) - pmap = Dict([pa => p[i] for pa => i in pidxmap]) + stmap = Dict([st => u[1][i] for (st, i) in stidxmap]) + pmap = Dict([pa => p[i] for (pa, i) in pidxmap]) eq = Symbolics.substitute(eq, merge(stmap, pmap)) csyms = [] # Find most nested calls, substitute those first. while !isempty(find_callable_syms!(csyms, eq)) for sym in csyms - t = arguments(sym)[1] x = operation(sym) + t = arguments(sym)[1] + prog = (tspan[2] - tspan[1])/(t - tspan[1]) # 1 / the % of the timespan elapsed if isparameter(x) eq = Symbolics.substitute(eq, Dict(x(t) => p[pidxmap[x(iv)]])) elseif isvariable(x) - eq = Symbolics.substitute(eq, Dict(x(t) => u(val)[stidxmap[x(iv)]])) + eq = Symbolics.substitute(eq, Dict(x(t) => u[Int(end ÷ prog)][stidxmap[x(iv)]])) end end empty!(csyms) diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 2d5535325a..6432c5ae02 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -1,80 +1,85 @@ -using BoundaryValueDiffEq, OrdinaryDiffEq +using BoundaryValueDiffEq, OrdinaryDiffEq, BoundaryValueDiffEqAscher using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D ### Test Collocation solvers on simple problems solvers = [MIRK4, RadauIIa5, LobattoIIIa3] +daesolvers = [Ascher2, Ascher4, Ascher6] -@parameters α=7.5 β=4.0 γ=8.0 δ=5.0 -@variables x(t)=1.0 y(t)=2.0 - -eqs = [D(x) ~ α * x - β * x * y, - D(y) ~ -γ * y + δ * x * y] - -u0map = [x => 1.0, y => 2.0] -parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0] -tspan = (0.0, 10.0) - -@mtkbuild lotkavolterra = ODESystem(eqs, t) -op = ODEProblem(lotkavolterra, u0map, tspan, parammap) -osol = solve(op, Vern9()) - -bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}( - lotkavolterra, u0map, tspan, parammap; eval_expression = true) - -for solver in solvers - sol = solve(bvp, solver(), dt = 0.01) - @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - @test sol.u[1] == [1.0, 2.0] -end - -# Test out of place -bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}( - lotkavolterra, u0map, tspan, parammap; eval_expression = true) - -for solver in solvers - sol = solve(bvp2, solver(), dt = 0.01) - @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - @test sol.u[1] == [1.0, 2.0] +let + @parameters α=7.5 β=4.0 γ=8.0 δ=5.0 + @variables x(t)=1.0 y(t)=2.0 + + eqs = [D(x) ~ α * x - β * x * y, + D(y) ~ -γ * y + δ * x * y] + + u0map = [x => 1.0, y => 2.0] + parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0] + tspan = (0.0, 10.0) + + @mtkbuild lotkavolterra = ODESystem(eqs, t) + op = ODEProblem(lotkavolterra, u0map, tspan, parammap) + osol = solve(op, Vern9()) + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}( + lotkavolterra, u0map, tspan, parammap; eval_expression = true) + + for solver in solvers + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [1.0, 2.0] + end + + # Test out of place + bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}( + lotkavolterra, u0map, tspan, parammap; eval_expression = true) + + for solver in solvers + sol = solve(bvp2, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [1.0, 2.0] + end end ### Testing on pendulum - -@parameters g=9.81 L=1.0 -@variables θ(t) = π / 2 - -eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] - -@mtkbuild pend = ODESystem(eqs, t) - -u0map = [θ => π / 2, D(θ) => π / 2] -parammap = [:L => 1.0, :g => 9.81] -tspan = (0.0, 6.0) - -op = ODEProblem(pend, u0map, tspan, parammap) -osol = solve(op, Vern9()) - -bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) -for solver in solvers - sol = solve(bvp, solver(), dt = 0.01) - @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - @test sol.u[1] == [π / 2, π / 2] -end - -# Test out-of-place -bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) - -for solver in solvers - sol = solve(bvp2, solver(), dt = 0.01) - @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - @test sol.u[1] == [π / 2, π / 2] +let + @parameters g=9.81 L=1.0 + @variables θ(t) = π / 2 + + eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] + + @mtkbuild pend = ODESystem(eqs, t) + + u0map = [θ => π / 2, D(θ) => π / 2] + parammap = [:L => 1.0, :g => 9.81] + tspan = (0.0, 6.0) + + op = ODEProblem(pend, u0map, tspan, parammap) + osol = solve(op, Vern9()) + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) + for solver in solvers + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [π / 2, π / 2] + end + + # Test out-of-place + bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) + + for solver in solvers + sol = solve(bvp2, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [π / 2, π / 2] + end end ################################################### -### TESTING ODESystem with Constraint Equations ### +### ODESystem with Constraint Equations, DAEs with constraints ### ################################################### -# Cartesian pendulum from the docs. Testing that initialization is satisfied. +# Cartesian pendulum from the docs. +# DAE IVP solved using BoundaryValueDiffEq solvers. let @parameters g @variables x(t) y(t) [state_priority = 10] λ(t) @@ -109,14 +114,74 @@ let end end -# Adding a midpoint boundary condition. +function test_solvers(solvers, prob, u0map, constraints, equations = []; dt = 0.01) + for solver in solvers + sol = solve(bvp, solver(); dt) + + for (k, v) in u0map + @test sol[k][1] == v + end + + for cons in constraints + @test sol[cons.rhs - cons.lhs] ≈ 0 + end + + for eq in equations + @test sol[eq] ≈ 0 + end + end +end + +# Simple ODESystem with BVP constraints. +let + @parameters α=7.5 β=4.0 γ=8.0 δ=5.0 + @variables x(..) y(t) + + eqs = [D(x) ~ α * x - β * x * y, + D(y) ~ -γ * y + δ * x * y] + + u0map = [y => 2.0] + parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0] + tspan = (0.0, 10.0) + guesses = [x => 1.0] + + @mtkbuild lotkavolterra = ODESystem(eqs, t) + op = ODEProblem(lotkavolterra, u0map, tspan, parammap, guesses = guesses) + + constraints = [x(6.) ~ 3] + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + test_solvers(solvers, bvp, u0map, constraints) + + # Testing that more complicated constraints give correct solutions. + constraints = [y(2.) + x(8.) ~ 12] + bvp = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + test_solvers(solvers, bvp, u0map, constraints) + + constraints = [α * β - x(6.) ~ 24] + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + test_solvers(solvers, bvp, u0map, constraints) + + # Testing that errors are properly thrown when malformed constraints are given. + @variables bad(..) + constraints = [x(1.) + bad(3.) ~ 10] + @test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + + constraints = [x(t) + y(t) ~ 3] + @test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + + @parameters bad2 + constraints = [bad2 + x(0.) ~ 3] + @test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) +end + +# Adding a midpoint boundary constraint. +# Solve using BVDAE solvers. let @parameters g @variables x(..) y(t) [state_priority = 10] λ(t) eqs = [D(D(x(t))) ~ λ * x(t) D(D(y)) ~ λ * y - g - x(t)^2 + y^2 ~ 1 - x(0.5) ~ 1] + x(t)^2 + y^2 ~ 1] @mtkbuild pend = ODESystem(eqs, t) tspan = (0.0, 1.5) @@ -124,37 +189,24 @@ let parammap = [g => 1] guesses = [λ => 1] - prob = ODEProblem(pend, u0map, tspan, pmap; guesses, check_length = false) - sol = solve(prob, Rodas5P()) + constraints = [x(0.5) ~ 1] + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) + test_solvers(daesolvers, bvp, u0map, constraints) - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses = guesses, check_length = false) - - for solver in solvers - sol = solve(bvp, solver(), dt = 0.01) - @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - conditions = getfield.(equations(pend)[3:end], :rhs) - @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 - @test sol.u[1] == [π / 2, π / 2] - end + bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) + test_solvers(daesolvers, bvp2, u0map, constraints, get_alg_eqs(pend)) - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses) - - for solver in solvers - sol = solve(bvp, solver(), dt = 0.01) - @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - conditions = getfield.(equations(pend)[3:end], :rhs) - @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 - end + # More complicated constraints. + u0map = [x(t) => 0.6] + guesses = [λ => 1, y(t) => 0.8] - bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) - for solver in solvers - sol = solve(bvp, solver(), dt = 0.01) - @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - conditions = getfield.(equations(pend)[3:end], :rhs) - @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 - end -end + constraints = [x(0.5) ~ 1, + x(0.3)^3 + y(0.6)^2 ~ 0.5] + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) + test_solvers(daesolvers, bvp, u0map, constraints, get_alg_eqs(pend)) -# Testing a more complicated case with multiple constraints. -let + constraints = [x(0.4) * g ~ y(0.2), + y(0.7) ~ 0.3] + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) + test_solvers(daesolvers, bvp, u0map, constraints, get_alg_eqs(pend)) end From 5d082ab05dc674b005ef1a724aa7273a712c2b66 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 11 Jan 2025 16:20:55 -0500 Subject: [PATCH 22/24] up --- src/systems/diffeqs/abstractodesystem.jl | 8 +++---- test/bvproblem.jl | 27 ++++++++++++------------ 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 3c45b1b2f8..99b34ec1a8 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -917,7 +917,7 @@ end function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], tspan = get_tspan(sys), parammap = DiffEqBase.NullParameters(); - constraints = nothing, guesses = nothing, + constraints = nothing, guesses = Dict(), version = nothing, tgrad = false, callback = nothing, check_length = true, @@ -929,7 +929,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") end - !isnothing(callbacks) && error("BVP solvers do not support callbacks.") + !isnothing(callback) && error("BVP solvers do not support callbacks.") iv = get_iv(sys) constraintsts = nothing @@ -964,7 +964,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] pidxmap = Dict([v => i for (i, v) in enumerate(ps)]) # Indices of states that have initial constraints. - u0i = has_alg_eqs(sys) ? collect(1:length(sts)) : [stidxmap[k] for k in keys(u0map)] + u0i = has_alg_eqs(sys) ? collect(1:length(sts)) : [stidxmap[k] for (k,v) in u0map] ni = length(u0i) bc = if !isnothing(constraints) @@ -994,7 +994,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] end end - return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...) + return BVProblem{iip}(f, bc, u0, tspan, p; kwargs...) end get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 6432c5ae02..e864433f3c 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -44,13 +44,14 @@ end ### Testing on pendulum let @parameters g=9.81 L=1.0 - @variables θ(t) = π / 2 + @variables θ(t) = π / 2 θ_t(t) - eqs = [D(D(θ)) ~ -(g / L) * sin(θ)] + eqs = [D(θ) ~ θ_t + D(θ_t) ~ -(g / L) * sin(θ)] @mtkbuild pend = ODESystem(eqs, t) - u0map = [θ => π / 2, D(θ) => π / 2] + u0map = [θ => π / 2, θ_t => π / 2] parammap = [:L => 1.0, :g => 9.81] tspan = (0.0, 6.0) @@ -74,9 +75,9 @@ let end end -################################################### -### ODESystem with Constraint Equations, DAEs with constraints ### -################################################### +################################################################## +### ODESystem with constraint equations, DAEs with constraints ### +################################################################## # Cartesian pendulum from the docs. # DAE IVP solved using BoundaryValueDiffEq solvers. @@ -90,19 +91,19 @@ let tspan = (0.0, 1.5) u0map = [x => 1, y => 0] - parammap = [g => 1] - guesses = [λ => 1] + pmap = [g => 1] + guess = [λ => 1] - prob = ODEProblem(pend, u0map, tspan, pmap; guesses) - sol = solve(prob, Rodas5P()) + prob = ODEProblem(pend, u0map, tspan, pmap; guesses = guess) + osol = solve(prob, Rodas5P()) - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses) + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses = guess) for solver in solvers - sol = solve(bvp, solver(), dt = 0.01) + sol = solve(bvp, solver(), dt = 0.001) @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) conditions = getfield.(equations(pend)[3:end], :rhs) - @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 + @test isapprox([sol[conditions][1]; sol[x][1] - 1; sol[y][1]], zeros(5), atol = 0.001) end bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) From b83e003babe1347439d05ea87a8b90995f1dcc82 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 14 Jan 2025 00:59:31 -0500 Subject: [PATCH 23/24] refactor the bc creation function --- src/systems/diffeqs/abstractodesystem.jl | 189 ++++++++++--------- test/bvproblem.jl | 230 ++++++++++++++++------- 2 files changed, 264 insertions(+), 155 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 99b34ec1a8..25347a17cd 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -931,68 +931,60 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] end !isnothing(callback) && error("BVP solvers do not support callbacks.") - iv = get_iv(sys) + has_alg_eqs(sys) && error("The BVProblem currently does not support ODESystems with algebraic equations.") # Remove this when the BVDAE solvers get updated, the codegen should work when it does. + constraintsts = nothing constraintps = nothing sts = unknowns(sys) ps = parameters(sys) - if !isnothing(constraints) + # Constraint validation + f_cons = if !isnothing(constraints) constraints isa Equation || constraints isa Vector{Equation} || error("Constraints must be specified as an equation or a vector of equations.") (length(constraints) + length(u0map) > length(sts)) && - error("The BVProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) cannot exceed the total number of states.") - - constraintsts = OrderedSet() - constraintps = OrderedSet() - - for eq in constraints - collect_vars!(constraintsts, constraintps, eq, iv) - validate_constraint_syms(eq, constraintsts, constraintps, Set(sts), Set(ps), iv) - empty!(constraintsts) - empty!(constraintps) - end + error("The BVProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) cannot exceed the total number of states.") end - f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; + # ODESystems without algebraic equations should use both fixed values + guesses + # for initialization. + _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses)) + f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, _u0map, parammap; t = tspan !== nothing ? tspan[1] : tspan, guesses, check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) stidxmap = Dict([v => i for (i, v) in enumerate(sts)]) - pidxmap = Dict([v => i for (i, v) in enumerate(ps)]) - - # Indices of states that have initial constraints. - u0i = has_alg_eqs(sys) ? collect(1:length(sts)) : [stidxmap[k] for (k,v) in u0map] - ni = length(u0i) - - bc = if !isnothing(constraints) - ne = length(constraints) - if iip - (residual,u,p,t) -> begin - residual[1:ni] .= u[1][u0i] .- u0[u0i] - residual[ni+1:ni+ne] .= map(constraints) do cons - sub_u_p_into_symeq(cons.rhs - cons.lhs, u, p, stidxmap, pidxmap, iv, tspan) - end - end - else - (u,p,t) -> begin - consresid = map(constraints) do cons - sub_u_p_into_symeq(cons.rhs-cons.lhs, u, p, stidxmap, pidxmap, iv, tspan) - end - resid = vcat(u[1][u0i] - u0[u0i], consresid) - end - end - else - if iip - (residual,u,p,t) -> begin - residual .= u[1] .- u0 - end - else - (u,p,t) -> (u[1] - u0) - end - end + u0_idxs = has_alg_eqs(sys) ? collect(1:length(sts)) : [stidxmap[k] for (k,v) in u0map] + + # bc = if !isnothing(constraints) && iip + # (residual,u,p,t) -> begin + # println(u(0.5)) + # residual[1:ni] .= u[1][u0i] .- u0[u0i] + # for (i, cons) in enumerate(constraints) + # residual[ni+i] = eval_symbolic_residual(cons, u, p, stidxmap, pidxmap, iv, tspan) + # end + # end + + # elseif !isnothing(constraints) && !iip + # (u,p,t) -> begin + # consresid = map(constraints) do cons + # eval_symbolic_residual(cons, u, p, stidxmap, pidxmap, iv, tspan) + # end + # resid = vcat(u[1][u0i] - u0[u0i], consresid) + # end + + # elseif iip + # (residual,u,p,t) -> begin + # println(u(0.5)) + # residual .= u[1] .- u0 + # end + + # else + # (u,p,t) -> (u[1] - u0) + # end + bc = process_constraints(sys, constraints, u0, u0_idxs, tspan, iip) return BVProblem{iip}(f, bc, u0, tspan, p; kwargs...) end @@ -1001,11 +993,10 @@ get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") # Validate that all the variables in the BVP constraints are well-formed states or parameters. function validate_constraint_syms(eq, constraintsts, constraintps, sts, ps, iv) - ModelingToolkit.check_variables(constraintsts) - ModelingToolkit.check_parameters(constraintps) - for var in constraintsts - if arguments(var) == iv + if length(arguments(var)) > 1 + error("Too many arguments for variable $var.") + elseif arguments(var) == iv var ∈ sts || error("Constraint equation $eq contains a variable $var that is not a variable of the ODESystem.") error("Constraint equation $eq contains a variable $var that does not have a specified argument. Such equations should be specified as algebraic equations to the ODESystem rather than a boundary constraints.") else @@ -1017,53 +1008,81 @@ function validate_constraint_syms(eq, constraintsts, constraintps, sts, ps, iv) if !iscall(var) var ∈ ps || error("Constraint equation $eq contains a parameter $var that is not a parameter of the ODESystem.") else + length(arguments(var)) > 1 && error("Too many arguments for parameter $var.") operation(var) ∈ ps || error("Constraint equations contain a parameter $var that is not a parameter of the ODESystem.") end end end -# Helper to substitute numeric values for u, p into the algebraic equations in the ODESystem. Used to construct the boundary condition function. -# Take a system with variables x,y, parameters g -# -# 1 + x(0) + y(0) → 1 + u[1][1] + u[1][2] -# x(0.5) → u(0.5)[1] -# x(0.5)*g(0.5) → u(0.5)[1]*p[1] -function sub_u_p_into_symeq(eq, u, p, stidxmap, pidxmap, iv, tspan) - eq = Symbolics.unwrap(eq) - - stmap = Dict([st => u[1][i] for (st, i) in stidxmap]) - pmap = Dict([pa => p[i] for (pa, i) in pidxmap]) - eq = Symbolics.substitute(eq, merge(stmap, pmap)) - - csyms = [] - # Find most nested calls, substitute those first. - while !isempty(find_callable_syms!(csyms, eq)) - for sym in csyms - x = operation(sym) - t = arguments(sym)[1] - prog = (tspan[2] - tspan[1])/(t - tspan[1]) # 1 / the % of the timespan elapsed - - if isparameter(x) - eq = Symbolics.substitute(eq, Dict(x(t) => p[pidxmap[x(iv)]])) - elseif isvariable(x) - eq = Symbolics.substitute(eq, Dict(x(t) => u[Int(end ÷ prog)][stidxmap[x(iv)]])) +""" + process_constraints(sys, constraints, u0, tspan, iip) + + Given an ODESystem with some constraints, generate the boundary condition function. +""" +function process_constraints(sys::ODESystem, constraints, u0, u0_idxs, tspan, iip) + + iv = get_iv(sys) + sts = get_unknowns(sys) + ps = get_ps(sys) + np = length(ps) + ns = length(sts) + + stidxmap = Dict([v => i for (i, v) in enumerate(sts)]) + pidxmap = Dict([v => i for (i, v) in enumerate(ps)]) + + @variables sol(..)[1:ns] p[1:np] + exprs = Any[] + + constraintsts = OrderedSet() + constraintps = OrderedSet() + + !isnothing(constraints) && for cons in constraints + collect_vars!(constraintsts, constraintps, cons, iv) + validate_constraint_syms(cons, constraintsts, constraintps, Set(sts), Set(ps), iv) + expr = cons.rhs - cons.lhs + + for st in constraintsts + x = operation(st) + t = arguments(st)[1] + idx = stidxmap[x(iv)] + + expr = Symbolics.substitute(expr, Dict(x(t) => sol(t)[idx])) + end + + for var in constraintps + if iscall(var) + x = operation(var) + t = arguments(var)[1] + idx = pidxmap[x] + + expr = Symbolics.substitute(expr, Dict(x(t) => p[idx])) + else + idx = pidxmap[var] + expr = Symbolics.substitute(expr, Dict(var => p[idx])) end end - empty!(csyms) + + empty!(constraintsts) + empty!(constraintps) + push!(exprs, expr) end - eq -end -function find_callable_syms!(csyms, ex) - ex = Symbolics.unwrap(ex) + init_cond_exprs = Any[] - if iscall(ex) - operation(ex) isa Symbolic && (arguments(ex)[1] isa Symbolic) && push!(csyms, ex) # only add leaf nodes - for arg in arguments(ex) - find_callable_syms!(csyms, arg) + for i in u0_idxs + expr = sol(tspan[1])[i] - u0[i] + push!(init_cond_exprs, expr) + end + + exprs = vcat(init_cond_exprs, exprs) + bcs = Symbolics.build_function(exprs, sol, p, expression = Val{false}) + if iip + return (resid, u, p, t) -> begin + bcs[2](resid, u, p) end + else + return (u, p, t) -> bcs[1](u, p) end - csyms end """ diff --git a/test/bvproblem.jl b/test/bvproblem.jl index e864433f3c..7e0d45b128 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -1,3 +1,5 @@ +### TODO: update when BoundaryValueDiffEqAscher is updated to use the normal boundary condition conventions + using BoundaryValueDiffEq, OrdinaryDiffEq, BoundaryValueDiffEqAscher using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D @@ -81,51 +83,140 @@ end # Cartesian pendulum from the docs. # DAE IVP solved using BoundaryValueDiffEq solvers. +# let +# @parameters g +# @variables x(t) y(t) [state_priority = 10] λ(t) +# eqs = [D(D(x)) ~ λ * x +# D(D(y)) ~ λ * y - g +# x^2 + y^2 ~ 1] +# @mtkbuild pend = ODESystem(eqs, t) +# +# tspan = (0.0, 1.5) +# u0map = [x => 1, y => 0] +# pmap = [g => 1] +# guess = [λ => 1] +# +# prob = ODEProblem(pend, u0map, tspan, pmap; guesses = guess) +# osol = solve(prob, Rodas5P()) +# +# zeta = [0., 0., 0., 0., 0.] +# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses = guess) +# +# for solver in solvers +# sol = solve(bvp, solver(zeta), dt = 0.001) +# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) +# conditions = getfield.(equations(pend)[3:end], :rhs) +# @test isapprox([sol[conditions][1]; sol[x][1] - 1; sol[y][1]], zeros(5), atol = 0.001) +# end +# +# bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) +# for solver in solvers +# sol = solve(bvp, solver(zeta), dt = 0.01) +# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) +# conditions = getfield.(equations(pend)[3:end], :rhs) +# @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 +# end +# end + +# Test generation of boundary condition function. let - @parameters g - @variables x(t) y(t) [state_priority = 10] λ(t) - eqs = [D(D(x)) ~ λ * x - D(D(y)) ~ λ * y - g - x^2 + y^2 ~ 1] - @mtkbuild pend = ODESystem(eqs, t) - - tspan = (0.0, 1.5) - u0map = [x => 1, y => 0] - pmap = [g => 1] - guess = [λ => 1] - - prob = ODEProblem(pend, u0map, tspan, pmap; guesses = guess) - osol = solve(prob, Rodas5P()) - - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses = guess) + @parameters α=7.5 β=4.0 γ=8.0 δ=5.0 + @variables x(..) y(t) + eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y, + D(y) ~ -γ * y + δ * x(t) * y] - for solver in solvers - sol = solve(bvp, solver(), dt = 0.001) - @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - conditions = getfield.(equations(pend)[3:end], :rhs) - @test isapprox([sol[conditions][1]; sol[x][1] - 1; sol[y][1]], zeros(5), atol = 0.001) + tspan = (0., 10.) + @mtkbuild lksys = ODESystem(eqs, t) + + function lotkavolterra!(du, u, p, t) + du[1] = p[1]*u[1] - p[2]*u[1]*u[2] + du[2] = -p[3]*u[2] + p[4]*u[1]*u[2] end - bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) - for solver in solvers - sol = solve(bvp, solver(), dt = 0.01) - @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) - conditions = getfield.(equations(pend)[3:end], :rhs) - @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 + function lotkavolterra(u, p, t) + [p[1]*u[1] - p[2]*u[1]*u[2], -p[3]*u[2] + p[4]*u[1]*u[2]] + end + # Compare the built bc function to the actual constructed one. + function bc!(resid, u, p, t) + resid[1] = u[1][1] - 1. + resid[2] = u[1][2] - 2. + nothing + end + function bc(u, p, t) + [u[1][1] - 1., u[1][2] - 2.] + end + + constraints = nothing + u0 = [1., 2.]; p = [7.5, 4., 8., 5.] + genbc_iip = ModelingToolkit.process_constraints(lksys, constraints, u0, [1, 2], tspan, true) + genbc_oop = ModelingToolkit.process_constraints(lksys, constraints, u0, [1, 2], tspan, false) + + bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, [1,2], tspan, p) + bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, [1,2], tspan, p) + + sol1 = solve(bvpi1, MIRK4(), dt = 0.01) + sol2 = solve(bvpi2, MIRK4(), dt = 0.01) + @test sol1 ≈ sol2 + + bvpo1 = BVProblem(lotkavolterra, bc, [1,2], tspan, p) + bvpo2 = BVProblem(lotkavolterra, genbc_oop, [1,2], tspan, p) + + sol1 = solve(bvpo1, MIRK4(), dt = 0.01) + sol2 = solve(bvpo2, MIRK4(), dt = 0.01) + @test sol1 ≈ sol2 + + # Test with a constraint. + constraints = [x(0.5) ~ 1.] + + function bc!(resid, u, p, t) + resid[1] = u[1][2] - 2. + resid[2] = u(0.5)[1] - 1. + end + function bc(u, p, t) + [u[1][2] - 2., u(0.5)[1] - 1.] end + + u0 = [1., 2.] + genbc_iip = ModelingToolkit.process_constraints(lksys, constraints, u0, [2], tspan, true) + genbc_oop = ModelingToolkit.process_constraints(lksys, constraints, u0, [2], tspan, false) + + bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, [1,2], tspan, p) + bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, [1,2], tspan, p) + + sol1 = solve(bvpi1, MIRK4(), dt = 0.01) + sol2 = solve(bvpi2, MIRK4(), dt = 0.01) + @test sol1 ≈ sol2 # don't get true equality here, not sure why + + bvpo1 = BVProblem(lotkavolterra, bc, [1,2], tspan, p) + bvpo2 = BVProblem(lotkavolterra, genbc_oop, [1,2], tspan, p) + + sol1 = solve(bvpo1, MIRK4(), dt = 0.01) + sol2 = solve(bvpo2, MIRK4(), dt = 0.01) + @test sol1 ≈ sol2 end function test_solvers(solvers, prob, u0map, constraints, equations = []; dt = 0.01) for solver in solvers - sol = solve(bvp, solver(); dt) + sol = solve(prob, solver(); dt) + @test successful_retcode(sol.retcode) + p = prob.p; t = sol.t; bc = prob.f.bc + ns = length(prob.u0) + + if isinplace(bvp.f) + resid = zeros(ns) + bc!(resid, sol, p, t) + @test isapprox(zeros(ns), resid) + else + @test isapprox(zeros(ns), bc(sol, p, t)) + end for (k, v) in u0map @test sol[k][1] == v end - for cons in constraints - @test sol[cons.rhs - cons.lhs] ≈ 0 - end + # for cons in constraints + # @test sol[cons.rhs - cons.lhs] ≈ 0 + # end for eq in equations @test sol[eq] ≈ 0 @@ -138,19 +229,18 @@ let @parameters α=7.5 β=4.0 γ=8.0 δ=5.0 @variables x(..) y(t) - eqs = [D(x) ~ α * x - β * x * y, - D(y) ~ -γ * y + δ * x * y] + eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y, + D(y) ~ -γ * y + δ * x(t) * y] u0map = [y => 2.0] parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0] tspan = (0.0, 10.0) - guesses = [x => 1.0] + guesses = [x(t) => 1.0] @mtkbuild lotkavolterra = ODESystem(eqs, t) - op = ODEProblem(lotkavolterra, u0map, tspan, parammap, guesses = guesses) constraints = [x(6.) ~ 3] - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; guesses, constraints) test_solvers(solvers, bvp, u0map, constraints) # Testing that more complicated constraints give correct solutions. @@ -177,37 +267,37 @@ end # Adding a midpoint boundary constraint. # Solve using BVDAE solvers. -let - @parameters g - @variables x(..) y(t) [state_priority = 10] λ(t) - eqs = [D(D(x(t))) ~ λ * x(t) - D(D(y)) ~ λ * y - g - x(t)^2 + y^2 ~ 1] - @mtkbuild pend = ODESystem(eqs, t) - - tspan = (0.0, 1.5) - u0map = [x(t) => 0.6, y => 0.8] - parammap = [g => 1] - guesses = [λ => 1] - - constraints = [x(0.5) ~ 1] - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) - test_solvers(daesolvers, bvp, u0map, constraints) - - bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) - test_solvers(daesolvers, bvp2, u0map, constraints, get_alg_eqs(pend)) - - # More complicated constraints. - u0map = [x(t) => 0.6] - guesses = [λ => 1, y(t) => 0.8] - - constraints = [x(0.5) ~ 1, - x(0.3)^3 + y(0.6)^2 ~ 0.5] - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) - test_solvers(daesolvers, bvp, u0map, constraints, get_alg_eqs(pend)) - - constraints = [x(0.4) * g ~ y(0.2), - y(0.7) ~ 0.3] - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) - test_solvers(daesolvers, bvp, u0map, constraints, get_alg_eqs(pend)) -end +# let +# @parameters g +# @variables x(..) y(t) [state_priority = 10] λ(t) +# eqs = [D(D(x(t))) ~ λ * x(t) +# D(D(y)) ~ λ * y - g +# x(t)^2 + y^2 ~ 1] +# @mtkbuild pend = ODESystem(eqs, t) +# +# tspan = (0.0, 1.5) +# u0map = [x(t) => 0.6, y => 0.8] +# parammap = [g => 1] +# guesses = [λ => 1] +# +# constraints = [x(0.5) ~ 1] +# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) +# test_solvers(daesolvers, bvp, u0map, constraints) +# +# bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) +# test_solvers(daesolvers, bvp2, u0map, constraints, get_alg_eqs(pend)) +# +# # More complicated constraints. +# u0map = [x(t) => 0.6] +# guesses = [λ => 1, y(t) => 0.8] +# +# constraints = [x(0.5) ~ 1, +# x(0.3)^3 + y(0.6)^2 ~ 0.5] +# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) +# test_solvers(daesolvers, bvp, u0map, constraints, get_alg_eqs(pend)) +# +# constraints = [x(0.4) * g ~ y(0.2), +# y(0.7) ~ 0.3] +# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) +# test_solvers(daesolvers, bvp, u0map, constraints, get_alg_eqs(pend)) +# end From db5eb66ea29533e51aefd6ea49ea04e0257f3201 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 14 Jan 2025 13:46:31 -0500 Subject: [PATCH 24/24] up --- src/systems/diffeqs/abstractodesystem.jl | 33 +----- test/bvproblem.jl | 144 ++++++++++++----------- 2 files changed, 76 insertions(+), 101 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 25347a17cd..0e8435942c 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -939,7 +939,7 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] ps = parameters(sys) # Constraint validation - f_cons = if !isnothing(constraints) + if !isnothing(constraints) constraints isa Equation || constraints isa Vector{Equation} || error("Constraints must be specified as an equation or a vector of equations.") @@ -958,32 +958,6 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [] stidxmap = Dict([v => i for (i, v) in enumerate(sts)]) u0_idxs = has_alg_eqs(sys) ? collect(1:length(sts)) : [stidxmap[k] for (k,v) in u0map] - # bc = if !isnothing(constraints) && iip - # (residual,u,p,t) -> begin - # println(u(0.5)) - # residual[1:ni] .= u[1][u0i] .- u0[u0i] - # for (i, cons) in enumerate(constraints) - # residual[ni+i] = eval_symbolic_residual(cons, u, p, stidxmap, pidxmap, iv, tspan) - # end - # end - - # elseif !isnothing(constraints) && !iip - # (u,p,t) -> begin - # consresid = map(constraints) do cons - # eval_symbolic_residual(cons, u, p, stidxmap, pidxmap, iv, tspan) - # end - # resid = vcat(u[1][u0i] - u0[u0i], consresid) - # end - - # elseif iip - # (residual,u,p,t) -> begin - # println(u(0.5)) - # residual .= u[1] .- u0 - # end - - # else - # (u,p,t) -> (u[1] - u0) - # end bc = process_constraints(sys, constraints, u0, u0_idxs, tspan, iip) return BVProblem{iip}(f, bc, u0, tspan, p; kwargs...) @@ -1075,11 +1049,10 @@ function process_constraints(sys::ODESystem, constraints, u0, u0_idxs, tspan, ii end exprs = vcat(init_cond_exprs, exprs) + @show exprs bcs = Symbolics.build_function(exprs, sol, p, expression = Val{false}) if iip - return (resid, u, p, t) -> begin - bcs[2](resid, u, p) - end + return (resid, u, p, t) -> bcs[2](resid, u, p) else return (u, p, t) -> bcs[1](u, p) end diff --git a/test/bvproblem.jl b/test/bvproblem.jl index 7e0d45b128..42762b8a3f 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -2,7 +2,9 @@ using BoundaryValueDiffEq, OrdinaryDiffEq, BoundaryValueDiffEqAscher using ModelingToolkit +using SciMLBase using ModelingToolkit: t_nounits as t, D_nounits as D +import ModelingToolkit: process_constraints ### Test Collocation solvers on simple problems solvers = [MIRK4, RadauIIa5, LobattoIIIa3] @@ -81,46 +83,9 @@ end ### ODESystem with constraint equations, DAEs with constraints ### ################################################################## -# Cartesian pendulum from the docs. -# DAE IVP solved using BoundaryValueDiffEq solvers. -# let -# @parameters g -# @variables x(t) y(t) [state_priority = 10] λ(t) -# eqs = [D(D(x)) ~ λ * x -# D(D(y)) ~ λ * y - g -# x^2 + y^2 ~ 1] -# @mtkbuild pend = ODESystem(eqs, t) -# -# tspan = (0.0, 1.5) -# u0map = [x => 1, y => 0] -# pmap = [g => 1] -# guess = [λ => 1] -# -# prob = ODEProblem(pend, u0map, tspan, pmap; guesses = guess) -# osol = solve(prob, Rodas5P()) -# -# zeta = [0., 0., 0., 0., 0.] -# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses = guess) -# -# for solver in solvers -# sol = solve(bvp, solver(zeta), dt = 0.001) -# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) -# conditions = getfield.(equations(pend)[3:end], :rhs) -# @test isapprox([sol[conditions][1]; sol[x][1] - 1; sol[y][1]], zeros(5), atol = 0.001) -# end -# -# bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) -# for solver in solvers -# sol = solve(bvp, solver(zeta), dt = 0.01) -# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) -# conditions = getfield.(equations(pend)[3:end], :rhs) -# @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 -# end -# end - # Test generation of boundary condition function. let - @parameters α=7.5 β=4.0 γ=8.0 δ=5.0 + @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 @variables x(..) y(t) eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y, D(y) ~ -γ * y + δ * x(t) * y] @@ -130,11 +95,11 @@ let function lotkavolterra!(du, u, p, t) du[1] = p[1]*u[1] - p[2]*u[1]*u[2] - du[2] = -p[3]*u[2] + p[4]*u[1]*u[2] + du[2] = -p[4]*u[2] + p[3]*u[1]*u[2] end function lotkavolterra(u, p, t) - [p[1]*u[1] - p[2]*u[1]*u[2], -p[3]*u[2] + p[4]*u[1]*u[2]] + [p[1]*u[1] - p[2]*u[1]*u[2], -p[4]*u[2] + p[3]*u[1]*u[2]] end # Compare the built bc function to the actual constructed one. function bc!(resid, u, p, t) @@ -146,23 +111,22 @@ let [u[1][1] - 1., u[1][2] - 2.] end - constraints = nothing - u0 = [1., 2.]; p = [7.5, 4., 8., 5.] - genbc_iip = ModelingToolkit.process_constraints(lksys, constraints, u0, [1, 2], tspan, true) - genbc_oop = ModelingToolkit.process_constraints(lksys, constraints, u0, [1, 2], tspan, false) + u0 = [1., 2.]; p = [1.5, 1., 3., 1.] + genbc_iip = ModelingToolkit.process_constraints(lksys, nothing, u0, [1, 2], tspan, true) + genbc_oop = ModelingToolkit.process_constraints(lksys, nothing, u0, [1, 2], tspan, false) - bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, [1,2], tspan, p) - bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, [1,2], tspan, p) + bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, [1.,2.], tspan, p) + bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, [1.,2.], tspan, p) - sol1 = solve(bvpi1, MIRK4(), dt = 0.01) - sol2 = solve(bvpi2, MIRK4(), dt = 0.01) + sol1 = solve(bvpi1, MIRK4(), dt = 0.05) + sol2 = solve(bvpi2, MIRK4(), dt = 0.05) @test sol1 ≈ sol2 bvpo1 = BVProblem(lotkavolterra, bc, [1,2], tspan, p) bvpo2 = BVProblem(lotkavolterra, genbc_oop, [1,2], tspan, p) - sol1 = solve(bvpo1, MIRK4(), dt = 0.01) - sol2 = solve(bvpo2, MIRK4(), dt = 0.01) + sol1 = solve(bvpo1, MIRK4(), dt = 0.05) + sol2 = solve(bvpo2, MIRK4(), dt = 0.05) @test sol1 ≈ sol2 # Test with a constraint. @@ -180,28 +144,28 @@ let genbc_iip = ModelingToolkit.process_constraints(lksys, constraints, u0, [2], tspan, true) genbc_oop = ModelingToolkit.process_constraints(lksys, constraints, u0, [2], tspan, false) - bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, [1,2], tspan, p) - bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, [1,2], tspan, p) - - sol1 = solve(bvpi1, MIRK4(), dt = 0.01) - sol2 = solve(bvpi2, MIRK4(), dt = 0.01) + bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, [1.,2.], tspan, p) + bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, [1.,2.], tspan, p) + bvpi3 = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, u0map, tspan, parammap; guesses, constraints) + + sol1 = solve(bvpi1, MIRK4(), dt = 0.05) + sol2 = solve(bvpi2, MIRK4(), dt = 0.05) @test sol1 ≈ sol2 # don't get true equality here, not sure why bvpo1 = BVProblem(lotkavolterra, bc, [1,2], tspan, p) bvpo2 = BVProblem(lotkavolterra, genbc_oop, [1,2], tspan, p) - sol1 = solve(bvpo1, MIRK4(), dt = 0.01) - sol2 = solve(bvpo2, MIRK4(), dt = 0.01) + sol1 = solve(bvpo1, MIRK4(), dt = 0.05) + sol2 = solve(bvpo2, MIRK4(), dt = 0.05) @test sol1 ≈ sol2 end -function test_solvers(solvers, prob, u0map, constraints, equations = []; dt = 0.01) +function test_solvers(solvers, prob, u0map, constraints, equations = []; dt = 0.05) for solver in solvers sol = solve(prob, solver(); dt) - @test successful_retcode(sol.retcode) + @test SciMLBase.successful_retcode(sol.retcode) p = prob.p; t = sol.t; bc = prob.f.bc ns = length(prob.u0) - if isinplace(bvp.f) resid = zeros(ns) bc!(resid, sol, p, t) @@ -226,45 +190,83 @@ end # Simple ODESystem with BVP constraints. let - @parameters α=7.5 β=4.0 γ=8.0 δ=5.0 + t = ModelingToolkit.t_nounits; D = ModelingToolkit.D_nounits + @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 @variables x(..) y(t) eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y, D(y) ~ -γ * y + δ * x(t) * y] - u0map = [y => 2.0] + u0map = [] parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0] tspan = (0.0, 10.0) - guesses = [x(t) => 1.0] + guesses = [x(t) => 1.0, y => 2.] @mtkbuild lotkavolterra = ODESystem(eqs, t) - constraints = [x(6.) ~ 3] + constraints = [x(6.) ~ 1.5] bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; guesses, constraints) test_solvers(solvers, bvp, u0map, constraints) # Testing that more complicated constraints give correct solutions. - constraints = [y(2.) + x(8.) ~ 12] - bvp = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + constraints = [y(2.) + x(8.) ~ 2.] + bvp = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(lotkavolterra, u0map, tspan, parammap; guesses, constraints) test_solvers(solvers, bvp, u0map, constraints) - constraints = [α * β - x(6.) ~ 24] - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + constraints = [α * β - x(6.) ~ 0.5] + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; guesses, constraints) test_solvers(solvers, bvp, u0map, constraints) # Testing that errors are properly thrown when malformed constraints are given. @variables bad(..) constraints = [x(1.) + bad(3.) ~ 10] - @test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + @test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; guesses, constraints) constraints = [x(t) + y(t) ~ 3] - @test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + @test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; guesses, constraints) @parameters bad2 constraints = [bad2 + x(0.) ~ 3] - @test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, constraints) + @test_throws Exception bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; guesses, constraints) end +# Cartesian pendulum from the docs. +# DAE IVP solved using BoundaryValueDiffEq solvers. +# let +# @parameters g +# @variables x(t) y(t) [state_priority = 10] λ(t) +# eqs = [D(D(x)) ~ λ * x +# D(D(y)) ~ λ * y - g +# x^2 + y^2 ~ 1] +# @mtkbuild pend = ODESystem(eqs, t) +# +# tspan = (0.0, 1.5) +# u0map = [x => 1, y => 0] +# pmap = [g => 1] +# guess = [λ => 1] +# +# prob = ODEProblem(pend, u0map, tspan, pmap; guesses = guess) +# osol = solve(prob, Rodas5P()) +# +# zeta = [0., 0., 0., 0., 0.] +# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses = guess) +# +# for solver in solvers +# sol = solve(bvp, solver(zeta), dt = 0.001) +# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) +# conditions = getfield.(equations(pend)[3:end], :rhs) +# @test isapprox([sol[conditions][1]; sol[x][1] - 1; sol[y][1]], zeros(5), atol = 0.001) +# end +# +# bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) +# for solver in solvers +# sol = solve(bvp, solver(zeta), dt = 0.01) +# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) +# conditions = getfield.(equations(pend)[3:end], :rhs) +# @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 +# end +# end + # Adding a midpoint boundary constraint. # Solve using BVDAE solvers. # let