From 96b310b8d738e320bdb8255c9756d28f46b52255 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 29 Jan 2024 15:37:30 +0530 Subject: [PATCH 1/3] feat!: deprecate ODAEProblem --- docs/src/basics/Composition.md | 6 +- .../modelingtoolkitize_index_reduction.md | 19 +----- docs/src/examples/tearing_parallelism.md | 2 +- docs/src/systems/ODESystem.md | 6 -- src/structural_transformation/codegen.jl | 47 +-------------- src/systems/diffeqs/odesystem.jl | 3 +- test/components.jl | 11 ++-- test/odaeproblem.jl | 59 ------------------- test/odesystem.jl | 4 +- test/runtests.jl | 1 - test/state_selection.jl | 4 +- test/structural_transformation/tearing.jl | 8 ++- test/symbolic_events.jl | 2 +- 13 files changed, 23 insertions(+), 149 deletions(-) delete mode 100644 test/odaeproblem.jl diff --git a/docs/src/basics/Composition.md b/docs/src/basics/Composition.md index dc7dce4225..d9806887ac 100644 --- a/docs/src/basics/Composition.md +++ b/docs/src/basics/Composition.md @@ -254,10 +254,10 @@ sol[reqn.R] ## Tearing Problem Construction -Some system types, specifically `ODESystem` and `NonlinearSystem`, can be further +Some system types (specifically `NonlinearSystem`) can be further reduced if `structural_simplify` has already been applied to them. This is done -by using the alternative problem constructors, `ODAEProblem` and `BlockNonlinearProblem` -respectively. In these cases, the constructor uses the knowledge of the +by using the alternative problem constructors (`BlockNonlinearProblem`). +In these cases, the constructor uses the knowledge of the strongly connected components calculated during the process of simplification as the basis for building pre-simplified nonlinear systems in the implicit solving. In summary: these problems are structurally modified, but could be diff --git a/docs/src/examples/modelingtoolkitize_index_reduction.md b/docs/src/examples/modelingtoolkitize_index_reduction.md index fea6cf19f3..5b7f3b329e 100644 --- a/docs/src/examples/modelingtoolkitize_index_reduction.md +++ b/docs/src/examples/modelingtoolkitize_index_reduction.md @@ -30,7 +30,7 @@ tspan = (0, 10.0) pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p) traced_sys = modelingtoolkitize(pendulum_prob) pendulum_sys = structural_simplify(dae_index_lowering(traced_sys)) -prob = ODAEProblem(pendulum_sys, [], tspan) +prob = ODEProblem(pendulum_sys, [], tspan) sol = solve(prob, Tsit5(), abstol = 1e-8, reltol = 1e-8) plot(sol, idxs = unknowns(traced_sys)) ``` @@ -162,20 +162,3 @@ variables which are symbolically eliminated, or any variable reordering done for enhanced parallelism/performance, still show up in the resulting plot and the plot is shown in the same order as the original numerical code. - -Note that we can even go a bit further. If we use the `ODAEProblem` -constructor, we can remove the algebraic equations from the unknowns of the -system and fully transform the index-3 DAE into an index-0 ODE which can -be solved via an explicit Runge-Kutta method: - -```@example indexred -traced_sys = modelingtoolkitize(pendulum_prob) -pendulum_sys = structural_simplify(dae_index_lowering(traced_sys)) -prob = ODAEProblem(pendulum_sys, Pair[], tspan) -sol = solve(prob, Tsit5(), abstol = 1e-8, reltol = 1e-8) -plot(sol, idxs = unknowns(traced_sys)) -``` - -And there you go: this has transformed the model from being too hard to -solve with implicit DAE solvers, to something that is easily solved with -explicit Runge-Kutta methods for non-stiff equations. diff --git a/docs/src/examples/tearing_parallelism.md b/docs/src/examples/tearing_parallelism.md index 1a79603528..3a6593909a 100644 --- a/docs/src/examples/tearing_parallelism.md +++ b/docs/src/examples/tearing_parallelism.md @@ -173,7 +173,7 @@ is that, your attempts to parallelize are neigh: performing parallelism after structural simplification greatly improves the problem that can be parallelized, so this is better than trying to do it by hand. -After performing this, you can construct the `ODEProblem`/`ODAEProblem` and set +After performing this, you can construct the `ODEProblem` and set `parallel_form` to use the exposed parallelism in multithreaded function constructions, but this showcases why `structural_simplify` is so important to that process. diff --git a/docs/src/systems/ODESystem.md b/docs/src/systems/ODESystem.md index 3f9a4dc45e..241b68b2b6 100644 --- a/docs/src/systems/ODESystem.md +++ b/docs/src/systems/ODESystem.md @@ -53,12 +53,6 @@ ODEProblem(sys::ModelingToolkit.AbstractODESystem, args...) SteadyStateProblem(sys::ModelingToolkit.AbstractODESystem, args...) ``` -## Torn Problem Constructors - -```@docs -ODAEProblem -``` - ## Expression Constructors ```@docs diff --git a/src/structural_transformation/codegen.jl b/src/structural_transformation/codegen.jl index d09e76d9d8..c306d21cd0 100644 --- a/src/structural_transformation/codegen.jl +++ b/src/structural_transformation/codegen.jl @@ -505,50 +505,7 @@ function build_observed_function(state, ts, var_eq_matching, var_sccs, expression ? ex : drop_expr(@RuntimeGeneratedFunction(ex)) end -""" - ODAEProblem{iip}(sys, u0map, tspan, parammap = DiffEqBase.NullParameters(); kw...) - -This constructor acts similar to the one for [`ODEProblem`](@ref) with the following changes: -`ODESystem`s can sometimes be further reduced if `structural_simplify` has -already been applied to them. -In these cases, the constructor uses the knowledge of the strongly connected -components calculated during the process of simplification as the basis for -building pre-simplified nonlinear systems in the implicit solving. - -In summary: these problems are structurally modified, but could be -more efficient and more stable. Note, the returned object is still of type -[`ODEProblem`](@ref). -""" struct ODAEProblem{iip} end -ODAEProblem(args...; kw...) = ODAEProblem{true}(args...; kw...) - -function ODAEProblem{iip}(sys, - u0map, - tspan, - parammap = DiffEqBase.NullParameters(); - callback = nothing, - use_union = true, - tofloat = true, - check = true, - kwargs...) where {iip} - eqs = equations(sys) - check && ModelingToolkit.check_operator_variables(eqs, Differential) - fun, dvs = build_torn_function(sys; kwargs...) - ps = parameters(sys) - defs = defaults(sys) - - defs = ModelingToolkit.mergedefaults(defs, parammap, ps) - defs = ModelingToolkit.mergedefaults(defs, u0map, dvs) - u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults = defs, tofloat = true) - p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults = defs, tofloat, use_union) - - cbs = process_events(sys; callback, kwargs...) - - kwargs = filter_kwargs(kwargs) - if cbs === nothing - ODEProblem{iip}(fun, u0, tspan, p; kwargs...) - else - ODEProblem{iip}(fun, u0, tspan, p; callback = cbs, kwargs...) - end -end +@deprecate ODAEProblem(args...; kw...) ODEProblem(args...; kw...) +@deprecate ODAEProblem{iip}(args...; kw...) where {iip} ODEProblem{iip}(args...; kw...) diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index b308c56895..880fdd1c6b 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -135,8 +135,7 @@ struct ODESystem <: AbstractODESystem """ discrete_subsystems::Any """ - A list of actual unknowns needed to be solved by solvers. Only - used for ODAEProblem. + A list of actual unknowns needed to be solved by solvers. """ solved_unknowns::Union{Nothing, Vector{Any}} """ diff --git a/test/components.jl b/test/components.jl index 7d19817e68..bca1cd9731 100644 --- a/test/components.jl +++ b/test/components.jl @@ -52,7 +52,7 @@ u0 = [capacitor.v => 0.0 prob = ODEProblem(sys, u0, (0, 10.0)) sol = solve(prob, Rodas4()) check_rc_sol(sol) -prob = ODAEProblem(sys, u0, (0, 10.0)) +prob = ODEProblem(sys, u0, (0, 10.0)) sol = solve(prob, Rodas4()) check_rc_sol(sol) @@ -80,7 +80,7 @@ let params = [param_r1 => 1.0, param_c1 => 1.0] tspan = (0.0, 10.0) - prob = ODAEProblem(sys, u0, tspan, params) + prob = ODEProblem(sys, u0, tspan, params) @test solve(prob, Tsit5()).retcode == ReturnCode.Success end @@ -97,7 +97,7 @@ let @named rc_model2 = compose(_rc_model2, [resistor, resistor2, capacitor, source, ground]) sys2 = structural_simplify(rc_model2) - prob2 = ODAEProblem(sys2, u0, (0, 10.0)) + prob2 = ODEProblem(sys2, u0, (0, 10.0)) sol2 = solve(prob2, Tsit5()) @test sol2[source.p.i] ≈ sol2[rc_model2.source.p.i] ≈ -sol2[capacitor.i] end @@ -138,7 +138,7 @@ sol_inner_outer = solve(prob, Rodas4()) u0 = [ capacitor.v => 0.0, ] -prob = ODAEProblem(sys, u0, (0, 10.0)) +prob = ODEProblem(sys, u0, (0, 10.0)) sol = solve(prob, Tsit5()) @test sol[resistor.p.i] == sol[capacitor.p.i] @@ -155,7 +155,6 @@ sys = structural_simplify(ll_model) @test length(equations(sys)) == 2 u0 = unknowns(sys) .=> 0 @test_nowarn ODEProblem(sys, u0, (0, 10.0)) -@test_nowarn ODAEProblem(sys, u0, (0, 10.0)) prob = DAEProblem(sys, Differential(t).(unknowns(sys)) .=> 0, u0, (0, 0.5)) sol = solve(prob, DFBDF()) @test sol.retcode == SciMLBase.ReturnCode.Success @@ -294,5 +293,5 @@ rc_eqs = [connect(capacitor.n, resistor.p) @named rc_model = compose(_rc_model, [resistor, capacitor, ground]) sys = structural_simplify(rc_model) -prob = ODAEProblem(sys, u0, (0, 10.0)) +prob = ODEProblem(sys, u0, (0, 10.0)) sol = solve(prob, Tsit5()) diff --git a/test/odaeproblem.jl b/test/odaeproblem.jl deleted file mode 100644 index c2de571cb0..0000000000 --- a/test/odaeproblem.jl +++ /dev/null @@ -1,59 +0,0 @@ -using ModelingToolkit, ModelingToolkitStandardLibrary, Test -using OrdinaryDiffEq -using ModelingToolkitStandardLibrary.Electrical -using ModelingToolkitStandardLibrary.Blocks - -function Segment(; name) - @named R = Resistor(; R = 1) - @named r = Resistor(; R = 1) - @named C = Capacitor(; C = 1) - - @named p1 = Pin() # top-left - @named p2 = Pin() # top-right - @named n = Pin() # bottom - - eqs = [connect(p1, R.p) - connect(R.n, p2, r.p) - connect(r.n, C.p) - connect(C.n, n)] - - return ODESystem(eqs, t, [], []; - name = name, - systems = [r, R, C, n, p1, p2]) -end - -function Strip(; name) - num_segments = 10 - # construct `num_segments` segments - segments = [Segment(; name = Symbol(:St_, seg)) - for seg in 1:num_segments] - - @named p1 = Pin() # top-left - @named p2 = Pin() # top-right - @named n = Pin() # bottom - - eqs = [connect(p1, segments[1].p1) - connect(p2, segments[end].p2) - [connect(n, seg.n) for seg in segments]... - [connect(segments[i].p2, segments[i + 1].p1) for i in 1:(num_segments - 1)]...] - - return ODESystem(eqs, t, [], []; name, - systems = [p1, p2, n, segments...]) -end - -@variables t -@named source = Voltage() -@named c = Constant(k = 0.01) - -@named ground = Ground() -@named strip = Strip() - -rc_eqs = [connect(c.output, source.V) - connect(source.p, strip.p1, strip.p2) - connect(strip.n, source.n, ground.g)] - -@named rc_model = ODESystem(rc_eqs, t, systems = [strip, c, source, ground]) -sys = structural_simplify(rc_model) - -prob = ODAEProblem(sys, [], (0, 10)) -@test_nowarn solve(prob, Tsit5()) diff --git a/test/odesystem.jl b/test/odesystem.jl index 1320fa7634..842341e483 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -925,7 +925,7 @@ let der = Differential(t) @named sys4 = ODESystem([der(x) ~ -y; der(y) ~ 1 + pp * y + x], t) sys4s = structural_simplify(sys4) - prob = ODAEProblem(sys4s, [x => 1.0, D(x) => 1.0], (0, 1.0)) + prob = ODEProblem(sys4s, [x => 1.0, D(x) => 1.0], (0, 1.0)) @test string.(unknowns(prob.f.sys)) == ["x(t)", "y(t)"] @test string.(parameters(prob.f.sys)) == ["pp"] @test string.(independent_variables(prob.f.sys)) == ["t"] @@ -978,7 +978,7 @@ let der = Differential(t) @named sys4 = ODESystem([der(x) ~ -y; der(y) ~ 1 + pp * y + x], t) sys4s = structural_simplify(sys4) - prob = ODAEProblem(sys4s, [x => 1.0, D(x) => 1.0], (0, 1.0)) + prob = ODEProblem(sys4s, [x => 1.0, D(x) => 1.0], (0, 1.0)) @test !isnothing(prob.f.sys) end diff --git a/test/runtests.jl b/test/runtests.jl index e9b392dd4a..ee4561864b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -35,7 +35,6 @@ end @safetestset "Constraints Test" include("constraints.jl") @safetestset "Reduction Test" include("reduction.jl") @safetestset "Split Parameters Test" include("split_parameters.jl") - @safetestset "ODAEProblem Test" include("odaeproblem.jl") @safetestset "StaticArrays Test" include("static_arrays.jl") @safetestset "Components Test" include("components.jl") @safetestset "Model Parsing Test" include("model_parsing.jl") diff --git a/test/state_selection.jl b/test/state_selection.jl index f4fbbb3c0b..d9835c5633 100644 --- a/test/state_selection.jl +++ b/test/state_selection.jl @@ -127,7 +127,7 @@ let D(return_pipe.fluid_port_a.m) => 0.0, D(supply_pipe.fluid_port_a.m) => 0.0] prob1 = ODEProblem(sys, u0, (0.0, 10.0), []) - prob2 = ODAEProblem(sys, u0, (0.0, 10.0), []) + prob2 = ODEProblem(sys, u0, (0.0, 10.0), []) prob3 = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, u0, (0.0, 10.0), []) @test solve(prob1, FBDF()).retcode == ReturnCode.Success #@test solve(prob2, FBDF()).retcode == ReturnCode.Success @@ -194,7 +194,7 @@ let mo_3 => 2 Ek_3 => 3] prob1 = ODEProblem(sys, u0, (0.0, 0.1)) - prob2 = ODAEProblem(sys, u0, (0.0, 0.1)) + prob2 = ODEProblem(sys, u0, (0.0, 0.1)) @test solve(prob1, FBDF()).retcode == ReturnCode.Success @test_broken solve(prob2, FBDF()).retcode == ReturnCode.Success end diff --git a/test/structural_transformation/tearing.jl b/test/structural_transformation/tearing.jl index 818daf6c95..6d168371f3 100644 --- a/test/structural_transformation/tearing.jl +++ b/test/structural_transformation/tearing.jl @@ -155,7 +155,9 @@ newdaesys = structural_simplify(daesys) @test equations(newdaesys) == [D(x) ~ z; 0 ~ y + sin(z) - p * t] @test equations(tearing_substitution(newdaesys)) == [D(x) ~ z; 0 ~ x + sin(z) - p * t] @test isequal(unknowns(newdaesys), [x, z]) -prob = ODAEProblem(newdaesys, [x => 1.0], (0, 1.0), [p => 0.2]) +@test isequal(states(newdaesys), [x, z]) +@test_deprecated ODAEProblem(newdaesys, [x => 1.0], (0, 1.0), [p => 0.2]) +prob = ODEProblem(newdaesys, [x => 1.0], (0, 1.0), [p => 0.2]) du = [0.0]; u = [1.0]; pr = 0.2; @@ -166,7 +168,7 @@ prob.f(du, u, pr, tt) # test the initial guess is respected @named sys = ODESystem(eqs, t, defaults = Dict(z => Inf)) -infprob = ODAEProblem(structural_simplify(sys), [x => 1.0], (0, 1.0), [p => 0.2]) +infprob = ODEProblem(structural_simplify(sys), [x => 1.0], (0, 1.0), [p => 0.2]) @test_throws Any infprob.f(du, u, pr, tt) sol1 = solve(prob, Tsit5()) @@ -214,6 +216,6 @@ u0 = [mass.s => 0.0 sys = structural_simplify(ms_model) @test ModelingToolkit.get_jac(sys)[] === ModelingToolkit.EMPTY_JAC @test ModelingToolkit.get_tgrad(sys)[] === ModelingToolkit.EMPTY_TGRAD -prob_complex = ODAEProblem(sys, u0, (0, 1.0)) +prob_complex = ODEProblem(sys, u0, (0, 1.0)) sol = solve(prob_complex, Tsit5()) @test all(sol[mass.v] .== 1) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index d85f2118fd..07d13eb164 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -289,7 +289,7 @@ eq = [vs ~ sin(2pi * t) ev = [sin(20pi * t) ~ 0.0] => [vmeasured ~ v] @named sys = ODESystem(eq, continuous_events = ev) sys = structural_simplify(sys) -prob = ODAEProblem(sys, zeros(2), (0.0, 5.1)) +prob = ODEProblem(sys, zeros(2), (0.0, 5.1)) sol = solve(prob, Tsit5()) @test all(minimum((0:0.1:5) .- sol.t', dims = 2) .< 0.0001) # test that the solver stepped every 0.1s as dictated by event @test sol([0.25])[vmeasured][] == sol([0.23])[vmeasured][] # test the hold property From 509aea8108e5fdb9bd951f91499cef58a1e0cd5a Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 29 Jan 2024 21:22:03 +0530 Subject: [PATCH 2/3] TEMP: change version to CI runs --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 2d327d1d0b..2a4c33f0ce 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.0.0" +version = "8.76.0" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" From 77c129050cf4399df78a695a5c61829358e3b45c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 1 Feb 2024 12:15:16 +0530 Subject: [PATCH 3/3] test: fix tests --- test/components.jl | 2 +- test/state_selection.jl | 2 +- test/structural_transformation/tearing.jl | 20 ++++++++++---------- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/test/components.jl b/test/components.jl index bca1cd9731..fbb7fe9e70 100644 --- a/test/components.jl +++ b/test/components.jl @@ -98,7 +98,7 @@ let [resistor, resistor2, capacitor, source, ground]) sys2 = structural_simplify(rc_model2) prob2 = ODEProblem(sys2, u0, (0, 10.0)) - sol2 = solve(prob2, Tsit5()) + sol2 = solve(prob2, Rosenbrock23()) @test sol2[source.p.i] ≈ sol2[rc_model2.source.p.i] ≈ -sol2[capacitor.i] end diff --git a/test/state_selection.jl b/test/state_selection.jl index d9835c5633..c203e70ad0 100644 --- a/test/state_selection.jl +++ b/test/state_selection.jl @@ -196,7 +196,7 @@ let prob1 = ODEProblem(sys, u0, (0.0, 0.1)) prob2 = ODEProblem(sys, u0, (0.0, 0.1)) @test solve(prob1, FBDF()).retcode == ReturnCode.Success - @test_broken solve(prob2, FBDF()).retcode == ReturnCode.Success + @test solve(prob2, FBDF()).retcode == ReturnCode.Success end let diff --git a/test/structural_transformation/tearing.jl b/test/structural_transformation/tearing.jl index 6d168371f3..a79a99080c 100644 --- a/test/structural_transformation/tearing.jl +++ b/test/structural_transformation/tearing.jl @@ -156,35 +156,35 @@ newdaesys = structural_simplify(daesys) @test equations(tearing_substitution(newdaesys)) == [D(x) ~ z; 0 ~ x + sin(z) - p * t] @test isequal(unknowns(newdaesys), [x, z]) @test isequal(states(newdaesys), [x, z]) -@test_deprecated ODAEProblem(newdaesys, [x => 1.0], (0, 1.0), [p => 0.2]) -prob = ODEProblem(newdaesys, [x => 1.0], (0, 1.0), [p => 0.2]) -du = [0.0]; -u = [1.0]; +@test_deprecated ODAEProblem(newdaesys, [x => 1.0, z => -0.5π], (0, 1.0), [p => 0.2]) +prob = ODEProblem(newdaesys, [x => 1.0, z => -0.5π], (0, 1.0), [p => 0.2]) +du = [0.0, 0.0]; +u = [1.0, -0.5π]; pr = 0.2; tt = 0.1; @test_skip (@ballocated $(prob.f)($du, $u, $pr, $tt)) == 0 prob.f(du, u, pr, tt) -@test du≈[-asin(u[1] - pr * tt)] atol=1e-5 +@test du ≈ [u[2], u[1] + sin(u[2]) - pr * tt] atol=1e-5 # test the initial guess is respected @named sys = ODESystem(eqs, t, defaults = Dict(z => Inf)) infprob = ODEProblem(structural_simplify(sys), [x => 1.0], (0, 1.0), [p => 0.2]) -@test_throws Any infprob.f(du, u, pr, tt) +@test_throws Any infprob.f(du, infprob.u0, pr, tt) -sol1 = solve(prob, Tsit5()) +sol1 = solve(prob, RosShamp4(), reltol=8e-7) sol2 = solve(ODEProblem{false}((u, p, t) -> [-asin(u[1] - pr * t)], [1.0], (0, 1.0), 0.2), Tsit5(), tstops = sol1.t, adaptive = false) -@test Array(sol1)≈Array(sol2) atol=1e-5 +@test Array(sol1[x])≈Array(sol2[1, :]) atol=1e-5 @test sol1[x] == first.(sol1.u) @test sol1[y] == first.(sol1.u) -@test sin.(sol1[z]) .+ sol1[y]≈pr[1] * sol1.t atol=1e-5 +@test sin.(sol1[z]) .+ sol1[y]≈pr[1] * sol1.t atol=5e-5 @test sol1[sin(z) + y]≈sin.(sol1[z]) .+ sol1[y] rtol=1e-12 @test sol1[y, :] == sol1[x, :] -@test (@. sin(sol1[z, :]) + sol1[y, :])≈pr * sol1.t atol=1e-5 +@test (@. sin(sol1[z, :]) + sol1[y, :])≈pr * sol1.t atol=5e-5 # 1426 function Translational_Mass(; name, m = 1.0)