Merge pull request #2435 from AayushSabharwal/as/deprecate-odaeproblem

feat!: deprecate ODAEProblem

Author: ChrisRackauckas

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

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.

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.

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

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...)

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}}
     """

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,8 +97,8 @@ let
     @named rc_model2 = compose(_rc_model2,
         [resistor, resistor2, capacitor, source, ground])
     sys2 = structural_simplify(rc_model2)
-    prob2 = ODAEProblem(sys2, u0, (0, 10.0))
-    sol2 = solve(prob2, Tsit5())
+    prob2 = ODEProblem(sys2, u0, (0, 10.0))
+    sol2 = solve(prob2, Rosenbrock23())
     @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())

test/odaeproblem.jl

@@ -935,7 +935,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"]
@@ -988,7 +988,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
 

test/odesystem.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")

test/runtests.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,9 +194,9 @@ 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
+    @test solve(prob2, FBDF()).retcode == ReturnCode.Success
 end
 
 let

test/state_selection.jl

@@ -155,34 +155,36 @@ 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])
-du = [0.0];
-u = [1.0];
+@test isequal(states(newdaesys), [x, z])
+@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 = ODAEProblem(structural_simplify(sys), [x => 1.0], (0, 1.0), [p => 0.2])
-@test_throws Any infprob.f(du, u, pr, tt)
+infprob = ODEProblem(structural_simplify(sys), [x => 1.0], (0, 1.0), [p => 0.2])
+@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)
@@ -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)

test/structural_transformation/tearing.jl

@@ -292,7 +292,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