From 5ca99ebb99b62fa666010de5b370bca3fefa156d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 29 Dec 2023 15:59:40 -0500 Subject: [PATCH 01/47] Fix up initializesystem for hierarchical models MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The first version only worked on non-hierarchical models. Now it handles all models. Here's some example usage. Here's an example model: ```julia using ModelingToolkit, DifferentialEquations @parameters t D = Differential(t) @connector Port begin p(t) dm(t)=0, [connect = Flow] end @connector Flange begin dx(t)=0 f(t), [connect = Flow] end # Components ---- @mtkmodel Orifice begin @parameters begin Cₒ=2.7 Aₒ=0.00094 ρ₀=1000 p′=0 end @variables begin dm(t)=0 p₁(t)=p′ p₂(t)=p′ end @components begin port₁ = Port(p=p′) port₂ = Port(p=p′) end begin u = dm/(ρ₀*Aₒ) end @equations begin dm ~ +port₁.dm dm ~ -port₂.dm p₁ ~ port₁.p p₂ ~ port₂.p p₁ - p₂ ~ (1/2)*ρ₀*u^2*Cₒ end end @mtkmodel Volume begin @parameters begin A=0.1 ρ₀=1000 β=2e9 direction=+1 p′ x′ end @variables begin p(t)=p′ x(t)=x′ dm(t)=0 f(t)=p′ * A dx(t)=0 r(t), [guess = 1000] dr(t), [guess = 1000] end @components begin port = Port(p=p′) flange = Flange(f=-p′ * A * direction) end @equations begin D(x) ~ dx D(r) ~ dr p ~ +port.p dm ~ +port.dm # mass is entering f ~ -flange.f * direction # force is leaving dx ~ flange.dx * direction r ~ ρ₀*(1 + p/β) dm ~ (r*dx*A) + (dr*x*A) f ~ p * A end end @mtkmodel Mass begin @parameters begin m = 100 f′ end @variables begin f(t)=f′ x(t)=0 dx(t)=0 ẍ(t)=f′/m end @components begin flange = Flange(f=f′) end @equations begin D(x) ~ dx D(dx) ~ ẍ f ~ flange.f dx ~ flange.dx m*ẍ ~ f end end @mtkmodel Actuator begin @parameters begin p₁′ p₂′ end begin #constants x′=0.5 A=0.1 end @components begin port₁ = Port(p=p₁′) port₂ = Port(p=p₂′) vol₁ = Volume(p′=p₁′, x′=x′, direction=-1) vol₂ = Volume(p′=p₂′, x′=x′, direction=+1) mass = Mass(f′=(p₂′ - p₁′)*A) flange = Flange(f=0) end @equations begin connect(port₁, vol₁.port) connect(port₂, vol₂.port) connect(vol₁.flange, vol₂.flange, mass.flange, flange) end end @mtkmodel Source begin @parameters begin p′ end @components begin port = Port(p=p′) end @equations begin port.p ~ p′ end end @mtkmodel Damper begin @parameters begin c = 1000 end @components begin flange = Flange(f=0) end @equations begin flange.f ~ c*flange.dx end end @mtkmodel System begin @components begin res₁ = Orifice(p′=300e5) res₂ = Orifice(p′=0) act = Actuator(p₁′=300e5, p₂′=0) src = Source(p′=300e5) snk = Source(p′=0) dmp = Damper() end @equations begin connect(src.port, res₁.port₁) connect(res₁.port₂, act.port₁) connect(act.port₂, res₂.port₁) connect(res₂.port₂, snk.port) connect(dmp.flange, act.flange) end end @mtkbuild sys = System() ``` It's having troubles initializing, so what is the system it's trying to initialize? ```julia initsys = initializesystem(sys) ``` That gives it symbolically. We have 98 equations for 54 variables, shesh that's overloaded. ```julia julia> initprob = NonlinearProblem(initsys,[]) ERROR: ArgumentError: Equations (98), states (54), and initial conditions (54) are of different lengths. To allow a different number of equations than states use kwarg check_length=false. Stacktrace: [1] check_eqs_u0(eqs::Vector{…}, dvs::Vector{…}, u0::Vector{…}; check_length::Bool, kwargs::@Kwargs{}) @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\abstractsystem.jl:1871 [2] process_NonlinearProblem(constructor::Type, sys::NonlinearSystem, u0map::Vector{…}, parammap::SciMLBase.NullParameters; version::Nothing, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, kwargs::@Kwargs{…}) @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:339 [3] process_NonlinearProblem @ C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:323 [inlined] [4] (NonlinearProblem{…})(sys::NonlinearSystem, u0map::Vector{…}, parammap::SciMLBase.NullParameters; check_length::Bool, kwargs::@Kwargs{}) @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:368 [5] NonlinearProblem (repeats 2 times) @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:365 [inlined] [6] #NonlinearProblem#698 @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:362 [inlined] [7] NonlinearProblem(sys::NonlinearSystem, args::Vector{Any}) @ ModelingToolkit C:\Users\accou\.julia\packages\ModelingToolkit\Gpzyo\src\systems\nonlinear\nonlinearsystem.jl:361 [8] top-level scope @ REPL[1]:1 Some type information was truncated. Use `show(err)` to see complete types. ``` So we can't build a NonlinearProblem. But we can do this: ```julia initprob = NonlinearProblem(initsys,[], check_length=false, checkbounds = true) lsqprob = NonlinearLeastSquaresProblem(NonlinearFunction(initprob.f.f, resid_prototype = zeros(98)), initprob.u0, initprob.p) initsol = solve(lsqprob, GaussNewton(), reltol = 1e-12, abstol = 1e-12) retcode: Success u: 54-element Vector{Float64}: 0.5 1015.0 0.5 1000.0 0.0 -1.208682818218319e-28 ⋮ 5.048709793414476e-28 -3.52499105861307e-29 -1.5146129380243427e-28 -3.0e6 -30000.0 -3.0e6 ``` And now we can make our initial conditions match that: ```julia allinit = states(initsys) .=> initsol prob = ODEProblem(sys, allinit, (0,0.1)) ``` and solve: ```julia sol = solve(prob) tmp = [0.0, 0.0, -0.0, -0.0, -0.0, -0.0, 8.123585990009498e-54, 1.1599794698483246e-51, 1.5375323209568473e-26, -2.914685306392616e-26] ┌ Warning: Instability detected. Aborting └ @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\3Mujd\src\integrator_interface.jl:602 ``` Here is modified the solver so tmp is the `resid` vector inside of the initializer. You can see it worked because the initialization errors are all zero. This is as opposed to with the default initial conditions: ```julia sysguesses = [ModelingToolkit.getguess(st) for st in states(sys)] hasaguess = findall(!isnothing, sysguesses) sysguesses = Dict(states(sys)[hasaguess] .=> sysguesses[hasaguess]) prob = ODEProblem(sys, sysguesses, (0,0.1)) sol = solve(prob) ``` ```julia julia> sol = solve(prob) tmp = [-0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3.0e7, 0.0, 50.0, 50.0] nlsol.resid = [-3.0e7, 0.0, 50.0, 50.0] nlsol.retcode = SciMLBase.ReturnCode.MaxIters retcode: InitialFailure Interpolation: specialized 4rd order "free" stiffness-aware interpolation t: 1-element Vector{Float64}: 0.0 u: 1-element Vector{Vector{Float64}}: [0.5, 1000.0, 0.5, 1000.0, 0.0, 0.0, 0.0, 0.0, 1000.0, 1000.0] ``` --- src/systems/nonlinear/initializesystem.jl | 70 +++++++++++------------ 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 49acdf0127..07f3d3d135 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -3,51 +3,51 @@ $(TYPEDSIGNATURES) Generate `NonlinearSystem` which initializes an ODE problem from specified initial conditions of an `ODESystem`. """ -function initializesystem(sys::ODESystem; name = nameof(sys), kwargs...) - if has_parent(sys) && (parent = get_parent(sys); parent !== nothing) - sys = parent - end +function initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(), kwargs...) sts, eqs = unknowns(sys), equations(sys) - idxs_diff = isdiffeq.(eqs) idxs_alge = .!idxs_diff - - # Algebraic equations and initial guesses are unchanged - eqs_ics = similar(eqs) - u0 = Vector{Any}(undef, length(sts)) - - eqs_ics[idxs_alge] .= eqs[idxs_alge] - u0[idxs_alge] .= getmetadata.(unwrap.(sts[idxs_alge]), - Symbolics.VariableDefaultValue, - nothing) - - for idx in findall(idxs_diff) - st = sts[idx] - if !hasdefault(st) - error("Invalid setup: unknown $(st) has no default value or equation.") - end - - def = getdefault(st) - if def isa Equation - if !hasguess(st) - error("Invalid setup: unknown $(st) has an initial condition equation with no guess.") + num_alge = sum(idxs_alge) + + # Start the equations list with algebraic equations + eqs_ics = eqs[idxs_alge] + u0 = Vector{Pair}(undef, 0) + defs = ModelingToolkit.defaults(sys) + + full_states = [sts;getfield.((observed(sys)),:lhs)] + + # Refactor to ODESystem construction + # should be ModelingToolkit.guesses(sys) + sysguesses = [ModelingToolkit.getguess(st) for st in full_states] + hasaguess = findall(!isnothing, sysguesses) + sysguesses = todict(full_states[hasaguess] .=> sysguesses[hasaguess]) + guesses = merge(sysguesses, todict(guesses)) + + for st in full_states + if st ∈ keys(defs) + def = defs[st] + + if def isa Equation + st ∉ keys(guesses) && error("Invalid setup: unknown $(st) has an initial condition equation with no guess.") + push!(eqs_ics,def) + push!(u0,st => guesses[st]) + else + push!(eqs_ics,st ~ def) + push!(u0, st => def) end - guess = getguess(st) - eqs_ics[idx] = def - - u0[idx] = guess + elseif st ∈ keys(guesses) + push!(u0,st => guesses[st]) else - eqs_ics[idx] = st ~ def - - u0[idx] = def + error("Invalid setup: unknown $(st) has no default value or initial guess") end end pars = parameters(sys) - sys_nl = NonlinearSystem(eqs_ics, - sts, + + sys_nl = NonlinearSystem([eqs_ics; observed(sys)], + full_states, pars; - defaults = Dict(sts .=> u0), + defaults = merge(ModelingToolkit.defaults(sys),todict(u0)), name, kwargs...) From 4717aee650670c4140adf9046571eb7d01a03bc9 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 2 Jan 2024 14:55:14 -0500 Subject: [PATCH 02/47] Handle overdetermined systems gracefully when `fully_determined = false` --- src/bipartite_graph.jl | 2 +- src/systems/nonlinear/initializesystem.jl | 18 ++++++++++-------- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/src/bipartite_graph.jl b/src/bipartite_graph.jl index 2c12b52c11..bd0f8af6b3 100644 --- a/src/bipartite_graph.jl +++ b/src/bipartite_graph.jl @@ -423,7 +423,7 @@ return `false` may not be matched. """ function maximal_matching(g::BipartiteGraph, srcfilter = vsrc -> true, dstfilter = vdst -> true, ::Type{U} = Unassigned) where {U} - matching = Matching{U}(ndsts(g)) + matching = Matching{U}(max(nsrcs(g), ndsts(g))) foreach(Iterators.filter(srcfilter, 𝑠vertices(g))) do vsrc construct_augmenting_path!(matching, g, vsrc, dstfilter) end diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 07f3d3d135..3c9c14459b 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -14,7 +14,7 @@ function initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(), u0 = Vector{Pair}(undef, 0) defs = ModelingToolkit.defaults(sys) - full_states = [sts;getfield.((observed(sys)),:lhs)] + full_states = [sts; getfield.((observed(sys)), :lhs)] # Refactor to ODESystem construction # should be ModelingToolkit.guesses(sys) @@ -28,26 +28,28 @@ function initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(), def = defs[st] if def isa Equation - st ∉ keys(guesses) && error("Invalid setup: unknown $(st) has an initial condition equation with no guess.") - push!(eqs_ics,def) - push!(u0,st => guesses[st]) + st ∉ keys(guesses) && + error("Invalid setup: unknown $(st) has an initial condition equation with no guess.") + push!(eqs_ics, def) + push!(u0, st => guesses[st]) else - push!(eqs_ics,st ~ def) + push!(eqs_ics, st ~ def) push!(u0, st => def) end elseif st ∈ keys(guesses) - push!(u0,st => guesses[st]) + push!(u0, st => guesses[st]) else error("Invalid setup: unknown $(st) has no default value or initial guess") end end pars = parameters(sys) + nleqs = [eqs_ics; observed(sys)] - sys_nl = NonlinearSystem([eqs_ics; observed(sys)], + sys_nl = NonlinearSystem(nleqs, full_states, pars; - defaults = merge(ModelingToolkit.defaults(sys),todict(u0)), + defaults = merge(ModelingToolkit.defaults(sys), todict(u0)), name, kwargs...) From 6d6d3fa104a63a1a6c8b859a7bbeaaa38c0e9b9b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 24 Feb 2024 11:09:54 -0500 Subject: [PATCH 03/47] Setup guesses(sys) and passing override dictionaries --- src/systems/abstractsystem.jl | 5 +++++ src/systems/diffeqs/odesystem.jl | 20 +++++++++++++++++--- src/systems/nonlinear/initializesystem.jl | 8 +++----- 3 files changed, 25 insertions(+), 8 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index f95cb1fe6f..2d8e80bf27 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -551,6 +551,7 @@ for prop in [:eqs :var_to_name :ctrls :defaults + :guesses :observed :tgrad :jac @@ -933,6 +934,10 @@ function full_parameters(sys::AbstractSystem) vcat(parameters(sys), dependent_parameters(sys)) end +function guesses(sys::AbstractSystem) + get_guesses(sys) +end + # required in `src/connectors.jl:437` parameters(_) = [] diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index fe52c032e0..3d0e3ef059 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -88,6 +88,11 @@ struct ODESystem <: AbstractODESystem """ defaults::Dict """ + The guesses to use as the initial conditions for the + initialization system. + """ + guesses::Dict + """ Tearing result specifying how to solve the system. """ torn_matching::Union{Matching, Nothing} @@ -157,7 +162,7 @@ struct ODESystem <: AbstractODESystem parent::Any function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, - jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, + jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, torn_matching, connector_type, preface, cevents, devents, parameter_dependencies, metadata = nothing, gui_metadata = nothing, tearing_state = nothing, @@ -175,7 +180,7 @@ struct ODESystem <: AbstractODESystem check_units(u, deqs) end new(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac, - ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching, + ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, torn_matching, connector_type, preface, cevents, devents, parameter_dependencies, metadata, gui_metadata, tearing_state, substitutions, complete, index_cache, discrete_subsystems, solved_unknowns, split_idxs, parent) @@ -191,6 +196,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), + guesses = Dict(), connector_type = nothing, preface = nothing, continuous_events = nothing, @@ -217,6 +223,13 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; var_to_name = Dict() process_variables!(var_to_name, defaults, dvs′) process_variables!(var_to_name, defaults, ps′) + + sysguesses = [ModelingToolkit.getguess(st) for st in dvs′] + hasaguess = findall(!isnothing, sysguesses) + var_guesses = dvs′[hasaguess] .=> sysguesses[hasaguess] + sysguesses = isempty(var_guesses) ? Dict() : todict(var_guesses) + guesses = merge(sysguesses, todict(guesses)) + isempty(observed) || collect_var_to_name!(var_to_name, (eq.lhs for eq in observed)) tgrad = RefValue(EMPTY_TGRAD) @@ -234,11 +247,12 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; parameter_dependencies, ps′) ODESystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)), deqs, iv′, dvs′, ps′, tspan, var_to_name, ctrl′, observed, tgrad, jac, - ctrl_jac, Wfact, Wfact_t, name, systems, defaults, nothing, + ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, nothing, connector_type, preface, cont_callbacks, disc_callbacks, parameter_dependencies, metadata, gui_metadata, checks = checks) end + function ODESystem(eqs, iv; kwargs...) eqs = collect(eqs) # NOTE: this assumes that the order of algebraic equations doesn't matter diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 3c9c14459b..45e1f9bbaa 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -12,16 +12,14 @@ function initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(), # Start the equations list with algebraic equations eqs_ics = eqs[idxs_alge] u0 = Vector{Pair}(undef, 0) - defs = ModelingToolkit.defaults(sys) + defs = defaults(sys) full_states = [sts; getfield.((observed(sys)), :lhs)] # Refactor to ODESystem construction # should be ModelingToolkit.guesses(sys) - sysguesses = [ModelingToolkit.getguess(st) for st in full_states] - hasaguess = findall(!isnothing, sysguesses) - sysguesses = todict(full_states[hasaguess] .=> sysguesses[hasaguess]) - guesses = merge(sysguesses, todict(guesses)) + + guesses = merge(get_guesses(sys), todict(guesses)) for st in full_states if st ∈ keys(defs) From 14403046f6393ca22f8ea4668487945a2e4b0aa8 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 24 Feb 2024 11:31:53 -0500 Subject: [PATCH 04/47] Add NonlinearLeastSquaresProblem building --- src/systems/nonlinear/nonlinearsystem.jl | 73 +++++++++++++++++++++++- 1 file changed, 72 insertions(+), 1 deletion(-) diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index b00efde2ab..f7df5a611d 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -289,6 +289,7 @@ function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(s NonlinearFunction{iip}(f, sys = sys, jac = _jac === nothing ? nothing : _jac, + resid_prototype = length(dvs) == length(equations(sys)) ? nothing : zeros(length(equations(sys))), jac_prototype = sparse ? similar(calculate_jacobian(sys, sparse = sparse), Float64) : nothing, @@ -333,12 +334,13 @@ function NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = unknowns(sys), end jp_expr = sparse ? :(similar($(get_jac(sys)[]), Float64)) : :nothing - + resid_expr = length(dvs) == length(equations(sys)) ? :nothing : :(zeros($(length(equations(sys))))) ex = quote f = $f jac = $_jac NonlinearFunction{$iip}(f, jac = jac, + resid_prototype = resid_expr, jac_prototype = $jp_expr) end !linenumbers ? Base.remove_linenums!(ex) : ex @@ -399,6 +401,35 @@ function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem, u0map, NonlinearProblem{iip}(f, u0, p, pt; filter_kwargs(kwargs)...) end +""" +```julia +DiffEqBase.NonlinearLeastSquaresProblem{iip}(sys::NonlinearSystem, u0map, + parammap = DiffEqBase.NullParameters(); + jac = false, sparse = false, + checkbounds = false, + linenumbers = true, parallel = SerialForm(), + kwargs...) where {iip} +``` + +Generates an NonlinearProblem from a NonlinearSystem and allows for automatically +symbolically calculating numerical enhancements. +""" +function DiffEqBase.NonlinearLeastSquaresProblem(sys::NonlinearSystem, args...; kwargs...) + NonlinearLeastSquaresProblem{true}(sys, args...; kwargs...) +end + +function DiffEqBase.NonlinearLeastSquaresProblem{iip}(sys::NonlinearSystem, u0map, + parammap = DiffEqBase.NullParameters(); + check_length = false, kwargs...) where {iip} + if !iscomplete(sys) + error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearLeastSquaresProblem`") + end + f, u0, p = process_NonlinearProblem(NonlinearFunction{iip}, sys, u0map, parammap; + check_length, kwargs...) + pt = something(get_metadata(sys), StandardNonlinearProblem()) + NonlinearLeastSquaresProblem{iip}(f, u0, p; filter_kwargs(kwargs)...) +end + """ ```julia DiffEqBase.NonlinearProblemExpr{iip}(sys::NonlinearSystem, u0map, @@ -439,6 +470,46 @@ function NonlinearProblemExpr{iip}(sys::NonlinearSystem, u0map, !linenumbers ? Base.remove_linenums!(ex) : ex end +""" +```julia +DiffEqBase.NonlinearLeastSquaresProblemExpr{iip}(sys::NonlinearSystem, u0map, + parammap = DiffEqBase.NullParameters(); + jac = false, sparse = false, + checkbounds = false, + linenumbers = true, parallel = SerialForm(), + kwargs...) where {iip} +``` + +Generates a Julia expression for a NonlinearProblem from a +NonlinearSystem and allows for automatically symbolically calculating +numerical enhancements. +""" +struct NonlinearLeastSquaresProblemExpr{iip} end + +function NonlinearLeastSquaresProblemExpr(sys::NonlinearSystem, args...; kwargs...) + NonlinearLeastSquaresProblemExpr{true}(sys, args...; kwargs...) +end + +function NonlinearLeastSquaresProblemExpr{iip}(sys::NonlinearSystem, u0map, + parammap = DiffEqBase.NullParameters(); + check_length = false, + kwargs...) where {iip} + if !iscomplete(sys) + error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearProblemExpr`") + end + f, u0, p = process_NonlinearProblem(NonlinearFunctionExpr{iip}, sys, u0map, parammap; + check_length, kwargs...) + linenumbers = get(kwargs, :linenumbers, true) + + ex = quote + f = $f + u0 = $u0 + p = $p + NonlinearLeastSquaresProblem(f, u0, p; $(filter_kwargs(kwargs)...)) + end + !linenumbers ? Base.remove_linenums!(ex) : ex +end + function flatten(sys::NonlinearSystem, noeqs = false) systems = get_systems(sys) if isempty(systems) From ee77fb2779585f9c24d9eb61f5021f1e1e707921 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 24 Feb 2024 12:16:53 -0500 Subject: [PATCH 05/47] Build the initialization system generation into structural simplification --- src/ModelingToolkit.jl | 9 ++++----- src/systems/abstractsystem.jl | 1 + src/systems/diffeqs/odesystem.jl | 13 +++++++++---- src/systems/nonlinear/initializesystem.jl | 6 +++--- src/systems/systemstructure.jl | 18 ++++++++++++++++-- 5 files changed, 33 insertions(+), 14 deletions(-) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 5bc2532214..c0645f52b1 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -137,19 +137,18 @@ include("systems/model_parsing.jl") include("systems/connectors.jl") include("systems/callbacks.jl") +include("systems/nonlinear/nonlinearsystem.jl") include("systems/diffeqs/odesystem.jl") include("systems/diffeqs/sdesystem.jl") include("systems/diffeqs/abstractodesystem.jl") +include("systems/nonlinear/modelingtoolkitize.jl") +include("systems/nonlinear/initializesystem.jl") include("systems/diffeqs/first_order_transform.jl") include("systems/diffeqs/modelingtoolkitize.jl") include("systems/diffeqs/basic_transformations.jl") include("systems/jumps/jumpsystem.jl") -include("systems/nonlinear/nonlinearsystem.jl") -include("systems/nonlinear/modelingtoolkitize.jl") -include("systems/nonlinear/initializesystem.jl") - include("systems/optimization/constraints_system.jl") include("systems/optimization/optimizationsystem.jl") include("systems/optimization/modelingtoolkitize.jl") @@ -253,7 +252,7 @@ export toexpr, get_variables export simplify, substitute export build_function export modelingtoolkitize -export initializesystem +export initializesystem, generate_initializesystem export @variables, @parameters, @constants, @brownian export @named, @nonamespace, @namespace, extend, compose, complete diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 2d8e80bf27..353e0c32b6 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -572,6 +572,7 @@ for prop in [:eqs :connections :preface :torn_matching + :initializesystem :tearing_state :substitutions :metadata diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 3d0e3ef059..fb0968f7de 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -97,6 +97,10 @@ struct ODESystem <: AbstractODESystem """ torn_matching::Union{Matching, Nothing} """ + The system for performing the initialization. + """ + initializesystem::Union{Nothing, NonlinearSystem} + """ Type of the system. """ connector_type::Any @@ -163,7 +167,7 @@ struct ODESystem <: AbstractODESystem function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, - torn_matching, connector_type, preface, cevents, + torn_matching, initializesystem, connector_type, preface, cevents, devents, parameter_dependencies, metadata = nothing, gui_metadata = nothing, tearing_state = nothing, substitutions = nothing, complete = false, index_cache = nothing, @@ -181,7 +185,7 @@ struct ODESystem <: AbstractODESystem end new(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, torn_matching, - connector_type, preface, cevents, devents, parameter_dependencies, metadata, + initializesystem, connector_type, preface, cevents, devents, parameter_dependencies, metadata, gui_metadata, tearing_state, substitutions, complete, index_cache, discrete_subsystems, solved_unknowns, split_idxs, parent) end @@ -197,6 +201,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), guesses = Dict(), + initializesystem = nothing, connector_type = nothing, preface = nothing, continuous_events = nothing, @@ -247,7 +252,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; parameter_dependencies, ps′) ODESystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)), deqs, iv′, dvs′, ps′, tspan, var_to_name, ctrl′, observed, tgrad, jac, - ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, nothing, + ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, nothing, initializesystem, connector_type, preface, cont_callbacks, disc_callbacks, parameter_dependencies, metadata, gui_metadata, checks = checks) end @@ -548,4 +553,4 @@ function add_accumulations(sys::ODESystem, vars::Vector{<:Pair}) @set! sys.eqs = [eqs; Equation[D(a) ~ v[2] for (a, v) in zip(avars, vars)]] @set! sys.unknowns = [get_unknowns(sys); avars] @set! sys.defaults = merge(get_defaults(sys), Dict(a => 0.0 for a in avars)) -end +end \ No newline at end of file diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 45e1f9bbaa..e3dc15e3fd 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -3,7 +3,7 @@ $(TYPEDSIGNATURES) Generate `NonlinearSystem` which initializes an ODE problem from specified initial conditions of an `ODESystem`. """ -function initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(), kwargs...) +function generate_initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(), check_defguess = false, kwargs...) sts, eqs = unknowns(sys), equations(sys) idxs_diff = isdiffeq.(eqs) idxs_alge = .!idxs_diff @@ -26,7 +26,7 @@ function initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(), def = defs[st] if def isa Equation - st ∉ keys(guesses) && + st ∉ keys(guesses) && check_defguess && error("Invalid setup: unknown $(st) has an initial condition equation with no guess.") push!(eqs_ics, def) push!(u0, st => guesses[st]) @@ -36,7 +36,7 @@ function initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(), end elseif st ∈ keys(guesses) push!(u0, st => guesses[st]) - else + elseif check_defguess error("Invalid setup: unknown $(st) has no default value or initial guess") end end diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 064ea1314d..08b80989c5 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -568,7 +568,7 @@ function merge_io(io, inputs) end function structural_simplify!(state::TearingState, io = nothing; simplify = false, - check_consistency = true, fully_determined = true, + check_consistency = true, fully_determined = true, warn_initialize_determined = true, kwargs...) if state.sys isa ODESystem ci = ModelingToolkit.ClockInference(state) @@ -615,7 +615,7 @@ function structural_simplify!(state::TearingState, io = nothing; simplify = fals end function _structural_simplify!(state::TearingState, io; simplify = false, - check_consistency = true, fully_determined = true, + check_consistency = true, fully_determined = true, warn_initialize_determined = true, kwargs...) check_consistency &= fully_determined has_io = io !== nothing @@ -635,5 +635,19 @@ function _structural_simplify!(state::TearingState, io; simplify = false, end fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)] @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns) + + if sys isa ODESystem + isys = ModelingToolkit.generate_initializesystem(sys) + !isempty(equations(isys)) && (isys = structural_simplify(isys; fully_determined = false)) + @set! sys.initializesystem = isys + neqs = length(equations(isys)) + nunknown = length(unknowns(isys)) + if warn_initialize_determined && neqs > nunknown + @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares" + end + if warn_initialize_determined && neqs < nunknown + @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares" + end + end ModelingToolkit.invalidate_cache!(sys), input_idxs end From 6f960db6813ae5e4ed07f44e8afedfa5ff32fd1d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 24 Feb 2024 12:44:48 -0500 Subject: [PATCH 06/47] InitializationProblem works --- src/systems/abstractsystem.jl | 3 + src/systems/diffeqs/abstractodesystem.jl | 55 ++++++ test/initializationsystem.jl | 240 +++++++++++++++++++++++ test/runtests.jl | 1 + 4 files changed, 299 insertions(+) create mode 100644 test/initializationsystem.jl diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 353e0c32b6..3a3f2aedf6 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -537,6 +537,9 @@ function complete(sys::AbstractSystem; split = true) if split && has_index_cache(sys) @set! sys.index_cache = IndexCache(sys) end + if isdefined(sys, :initializationsystem) + @set! sys.initializationsystem = complete(get_initializationsystem(sys); split) + end isdefined(sys, :complete) ? (@set! sys.complete = true) : sys end diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index b6e8c50563..430e9b60a3 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1442,3 +1442,58 @@ function flatten_equations(eqs) end end end + +struct InitializationProblem{iip, specialization} end + +""" +```julia +InitializationProblem{iip}(sys::AbstractODESystem, u0map, tspan, + parammap = DiffEqBase.NullParameters(); + version = nothing, tgrad = false, + jac = false, + checkbounds = false, sparse = false, + simplify = false, + linenumbers = true, parallel = SerialForm(), + kwargs...) where {iip} +``` + +Generates a NonlinearProblem or NonlinearLeastSquaresProblem from an ODESystem +which represents the initialization, i.e. the calculation of the consistent +initial conditions for the given DAE. +""" +function InitializationProblem(sys::AbstractODESystem, args...; kwargs...) + InitializationProblem{true}(sys, args...; kwargs...) +end + +function InitializationProblem(sys::AbstractODESystem, + u0map::StaticArray, + args...; + kwargs...) + InitializationProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) +end + +function InitializationProblem{true}(sys::AbstractODESystem, args...; kwargs...) + InitializationProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) +end + +function InitializationProblem{false}(sys::AbstractODESystem, args...; kwargs...) + InitializationProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) +end + +function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], + parammap = DiffEqBase.NullParameters(); + check_length = true, + kwargs...) where {iip, specialize} + if !iscomplete(sys) + error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`") + end + + isys = get_initializesystem(sys) + neqs = length(equations(isys)) + nunknown = length(unknowns(isys)) + if neqs == nunknown + NonlinearProblem(isys,u0map,parammap) + else + NonlinearLeastSquaresProblem(isys,u0map,parammap) + end +end \ No newline at end of file diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl new file mode 100644 index 0000000000..d9e3c071a4 --- /dev/null +++ b/test/initializationsystem.jl @@ -0,0 +1,240 @@ +using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve, Test +using ModelingToolkit: t_nounits as t, D_nounits as D + +@connector Port begin + p(t) + dm(t)=0, [connect = Flow] +end + +@connector Flange begin + dx(t)=0 + f(t), [connect = Flow] +end + +# Components ---- +@mtkmodel Orifice begin + @parameters begin + Cₒ=2.7 + Aₒ=0.00094 + ρ₀=1000 + p′=0 + end + @variables begin + dm(t)=0 + p₁(t)=p′ + p₂(t)=p′ + end + @components begin + port₁ = Port(p=p′) + port₂ = Port(p=p′) + end + begin + u = dm/(ρ₀*Aₒ) + end + @equations begin + dm ~ +port₁.dm + dm ~ -port₂.dm + p₁ ~ port₁.p + p₂ ~ port₂.p + + p₁ - p₂ ~ (1/2)*ρ₀*u^2*Cₒ + end +end + +@mtkmodel Volume begin + @parameters begin + A=0.1 + ρ₀=1000 + β=2e9 + direction=+1 + p′ + x′ + end + @variables begin + p(t)=p′ + x(t)=x′ + dm(t)=0 + f(t)=p′ * A + dx(t)=0 + r(t), [guess = 1000] + dr(t), [guess = 1000] + end + @components begin + port = Port(p=p′) + flange = Flange(f=-p′ * A * direction) + end + @equations begin + D(x) ~ dx + D(r) ~ dr + + p ~ +port.p + dm ~ +port.dm # mass is entering + f ~ -flange.f * direction # force is leaving + dx ~ flange.dx * direction + + r ~ ρ₀*(1 + p/β) + dm ~ (r*dx*A) + (dr*x*A) + f ~ p * A + end +end + +@mtkmodel Mass begin + @parameters begin + m = 100 + f′ + end + @variables begin + f(t)=f′ + x(t)=0 + dx(t)=0 + ẍ(t)=f′/m + end + @components begin + flange = Flange(f=f′) + end + @equations begin + D(x) ~ dx + D(dx) ~ ẍ + + f ~ flange.f + dx ~ flange.dx + + m*ẍ ~ f + end +end + +@mtkmodel Actuator begin + @parameters begin + p₁′ + p₂′ + end + begin #constants + x′=0.5 + A=0.1 + end + @components begin + port₁ = Port(p=p₁′) + port₂ = Port(p=p₂′) + vol₁ = Volume(p′=p₁′, x′=x′, direction=-1) + vol₂ = Volume(p′=p₂′, x′=x′, direction=+1) + mass = Mass(f′=(p₂′ - p₁′)*A) + flange = Flange(f=0) + end + @equations begin + connect(port₁, vol₁.port) + connect(port₂, vol₂.port) + connect(vol₁.flange, vol₂.flange, mass.flange, flange) + end +end + +@mtkmodel Source begin + @parameters begin + p′ + end + @components begin + port = Port(p=p′) + end + @equations begin + port.p ~ p′ + end +end + +@mtkmodel Damper begin + @parameters begin + c = 1000 + end + @components begin + flange = Flange(f=0) + end + @equations begin + flange.f ~ c*flange.dx + end +end + +@mtkmodel System begin + @components begin + res₁ = Orifice(p′=300e5) + res₂ = Orifice(p′=0) + act = Actuator(p₁′=300e5, p₂′=0) + src = Source(p′=300e5) + snk = Source(p′=0) + dmp = Damper() + end + @equations begin + connect(src.port, res₁.port₁) + connect(res₁.port₂, act.port₁) + connect(act.port₂, res₂.port₁) + connect(res₂.port₂, snk.port) + connect(dmp.flange, act.flange) + end +end + +@mtkbuild sys = System() +initprob = ModelingToolkit.InitializationProblem(sys) +@test initprob isa NonlinearLeastSquaresProblem +@test length(initprob.u0) == 2 +initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12) +@test SciMLBase.successful_retcode(initsol) + +@connector Flange begin + dx(t), [guess = 0] + f(t), [guess = 0, connect=Flow] +end + +@mtkmodel Mass begin + @parameters begin + m = 100 + end + @variables begin + dx(t) + f(t), [guess = 0] + end + @components begin + flange = Flange() + end + @equations begin + # connectors + flange.dx ~ dx + flange.f ~ -f + + # physics + f ~ m*D(dx) + end +end + +@mtkmodel Damper begin + @parameters begin + d = 1 + end + @variables begin + dx(t), [guess = 0] + f(t), [guess = 0] + end + @components begin + flange = Flange() + end + @equations begin + # connectors + flange.dx ~ dx + flange.f ~ -f + + # physics + f ~ d*dx + end +end + +@mtkmodel MassDamperSystem begin + @components begin + mass = Mass(;dx=100,m=10) + damper = Damper(;d=1) + end + @equations begin + connect(mass.flange, damper.flange) + end +end + +@mtkbuild sys = MassDamperSystem() +initprob = ModelingToolkit.InitializationProblem(sys) +@test initprob isa NonlinearProblem +initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12) +@test SciMLBase.successful_retcode(initsol) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 96bdbbf06e..605825bc0b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,6 +34,7 @@ end @safetestset "SteadyStateSystem Test" include("steadystatesystems.jl") @safetestset "SDESystem Test" include("sdesystem.jl") @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl") + @safetestset "InitializationSystem Test" include("initializationsystem.jl") @safetestset "PDE Construction Test" include("pde.jl") @safetestset "JumpSystem Test" include("jumpsystem.jl") @safetestset "Constraints Test" include("constraints.jl") From 0f3ba28d58a96f178547a4dc6a74baf83268728c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 24 Feb 2024 13:04:22 -0500 Subject: [PATCH 07/47] Test the initialization problem --- test/initializationsystem.jl | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index d9e3c071a4..03caa0ddeb 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -176,6 +176,12 @@ initprob = ModelingToolkit.InitializationProblem(sys) initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12) @test SciMLBase.successful_retcode(initsol) +allinit = unknowns(sys) .=> initsol[unknowns(sys)] +prob = ODEProblem(sys, allinit, (0,0.1)) +sol = solve(prob, Rodas5P()) +# If initialized incorrectly, then it would be InitialFailure +@test sol.retcode == SciMLBase.ReturnCode.Unstable + @connector Flange begin dx(t), [guess = 0] f(t), [guess = 0, connect=Flow] @@ -237,4 +243,10 @@ end initprob = ModelingToolkit.InitializationProblem(sys) @test initprob isa NonlinearProblem initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12) -@test SciMLBase.successful_retcode(initsol) \ No newline at end of file +@test SciMLBase.successful_retcode(initsol) + +allinit = unknowns(sys) .=> initsol[unknowns(sys)] +prob = ODEProblem(sys, allinit, (0,0.1)) +sol = solve(prob, Rodas5P()) +# If initialized incorrectly, then it would be InitialFailure +@test sol.retcode == SciMLBase.ReturnCode.Success \ No newline at end of file From aafddc906c1bb8dfd435a49447a07187543ca981 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 24 Feb 2024 13:08:07 -0500 Subject: [PATCH 08/47] format --- src/systems/diffeqs/abstractodesystem.jl | 6 +- src/systems/diffeqs/odesystem.jl | 3 +- src/systems/nonlinear/initializesystem.jl | 5 +- src/systems/nonlinear/nonlinearsystem.jl | 6 +- src/systems/systemstructure.jl | 3 +- test/initializationsystem.jl | 116 +++++++++++----------- 6 files changed, 71 insertions(+), 68 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 430e9b60a3..42667101fe 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1492,8 +1492,8 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0map = neqs = length(equations(isys)) nunknown = length(unknowns(isys)) if neqs == nunknown - NonlinearProblem(isys,u0map,parammap) + NonlinearProblem(isys, u0map, parammap) else - NonlinearLeastSquaresProblem(isys,u0map,parammap) + NonlinearLeastSquaresProblem(isys, u0map, parammap) end -end \ No newline at end of file +end diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index fb0968f7de..7ee5f705f4 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -257,7 +257,6 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; metadata, gui_metadata, checks = checks) end - function ODESystem(eqs, iv; kwargs...) eqs = collect(eqs) # NOTE: this assumes that the order of algebraic equations doesn't matter @@ -553,4 +552,4 @@ function add_accumulations(sys::ODESystem, vars::Vector{<:Pair}) @set! sys.eqs = [eqs; Equation[D(a) ~ v[2] for (a, v) in zip(avars, vars)]] @set! sys.unknowns = [get_unknowns(sys); avars] @set! sys.defaults = merge(get_defaults(sys), Dict(a => 0.0 for a in avars)) -end \ No newline at end of file +end diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index e3dc15e3fd..80384bde4c 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -3,7 +3,8 @@ $(TYPEDSIGNATURES) Generate `NonlinearSystem` which initializes an ODE problem from specified initial conditions of an `ODESystem`. """ -function generate_initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(), check_defguess = false, kwargs...) +function generate_initializesystem(sys::ODESystem; name = nameof(sys), + guesses = Dict(), check_defguess = false, kwargs...) sts, eqs = unknowns(sys), equations(sys) idxs_diff = isdiffeq.(eqs) idxs_alge = .!idxs_diff @@ -18,7 +19,7 @@ function generate_initializesystem(sys::ODESystem; name = nameof(sys), guesses = # Refactor to ODESystem construction # should be ModelingToolkit.guesses(sys) - + guesses = merge(get_guesses(sys), todict(guesses)) for st in full_states diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index f7df5a611d..798ed10e8d 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -289,7 +289,8 @@ function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(s NonlinearFunction{iip}(f, sys = sys, jac = _jac === nothing ? nothing : _jac, - resid_prototype = length(dvs) == length(equations(sys)) ? nothing : zeros(length(equations(sys))), + resid_prototype = length(dvs) == length(equations(sys)) ? nothing : + zeros(length(equations(sys))), jac_prototype = sparse ? similar(calculate_jacobian(sys, sparse = sparse), Float64) : nothing, @@ -334,7 +335,8 @@ function NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = unknowns(sys), end jp_expr = sparse ? :(similar($(get_jac(sys)[]), Float64)) : :nothing - resid_expr = length(dvs) == length(equations(sys)) ? :nothing : :(zeros($(length(equations(sys))))) + resid_expr = length(dvs) == length(equations(sys)) ? :nothing : + :(zeros($(length(equations(sys))))) ex = quote f = $f jac = $_jac diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 08b80989c5..f18f5c5c89 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -638,7 +638,8 @@ function _structural_simplify!(state::TearingState, io; simplify = false, if sys isa ODESystem isys = ModelingToolkit.generate_initializesystem(sys) - !isempty(equations(isys)) && (isys = structural_simplify(isys; fully_determined = false)) + !isempty(equations(isys)) && + (isys = structural_simplify(isys; fully_determined = false)) @set! sys.initializesystem = isys neqs = length(equations(isys)) nunknown = length(unknowns(isys)) diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index 03caa0ddeb..a92698bc1a 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -3,33 +3,33 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @connector Port begin p(t) - dm(t)=0, [connect = Flow] + dm(t) = 0, [connect = Flow] end @connector Flange begin - dx(t)=0 + dx(t) = 0 f(t), [connect = Flow] end # Components ---- @mtkmodel Orifice begin @parameters begin - Cₒ=2.7 - Aₒ=0.00094 - ρ₀=1000 - p′=0 + Cₒ = 2.7 + Aₒ = 0.00094 + ρ₀ = 1000 + p′ = 0 end @variables begin - dm(t)=0 - p₁(t)=p′ - p₂(t)=p′ + dm(t) = 0 + p₁(t) = p′ + p₂(t) = p′ end @components begin - port₁ = Port(p=p′) - port₂ = Port(p=p′) + port₁ = Port(p = p′) + port₂ = Port(p = p′) end begin - u = dm/(ρ₀*Aₒ) + u = dm / (ρ₀ * Aₒ) end @equations begin dm ~ +port₁.dm @@ -37,31 +37,31 @@ end p₁ ~ port₁.p p₂ ~ port₂.p - p₁ - p₂ ~ (1/2)*ρ₀*u^2*Cₒ + p₁ - p₂ ~ (1 / 2) * ρ₀ * u^2 * Cₒ end end @mtkmodel Volume begin @parameters begin - A=0.1 - ρ₀=1000 - β=2e9 - direction=+1 + A = 0.1 + ρ₀ = 1000 + β = 2e9 + direction = +1 p′ x′ end @variables begin - p(t)=p′ - x(t)=x′ - dm(t)=0 - f(t)=p′ * A - dx(t)=0 + p(t) = p′ + x(t) = x′ + dm(t) = 0 + f(t) = p′ * A + dx(t) = 0 r(t), [guess = 1000] dr(t), [guess = 1000] end @components begin - port = Port(p=p′) - flange = Flange(f=-p′ * A * direction) + port = Port(p = p′) + flange = Flange(f = -p′ * A * direction) end @equations begin D(x) ~ dx @@ -72,8 +72,8 @@ end f ~ -flange.f * direction # force is leaving dx ~ flange.dx * direction - r ~ ρ₀*(1 + p/β) - dm ~ (r*dx*A) + (dr*x*A) + r ~ ρ₀ * (1 + p / β) + dm ~ (r * dx * A) + (dr * x * A) f ~ p * A end end @@ -84,13 +84,13 @@ end f′ end @variables begin - f(t)=f′ - x(t)=0 - dx(t)=0 - ẍ(t)=f′/m + f(t) = f′ + x(t) = 0 + dx(t) = 0 + ẍ(t) = f′ / m end @components begin - flange = Flange(f=f′) + flange = Flange(f = f′) end @equations begin D(x) ~ dx @@ -99,7 +99,7 @@ end f ~ flange.f dx ~ flange.dx - m*ẍ ~ f + m * ẍ ~ f end end @@ -109,16 +109,16 @@ end p₂′ end begin #constants - x′=0.5 - A=0.1 + x′ = 0.5 + A = 0.1 end @components begin - port₁ = Port(p=p₁′) - port₂ = Port(p=p₂′) - vol₁ = Volume(p′=p₁′, x′=x′, direction=-1) - vol₂ = Volume(p′=p₂′, x′=x′, direction=+1) - mass = Mass(f′=(p₂′ - p₁′)*A) - flange = Flange(f=0) + port₁ = Port(p = p₁′) + port₂ = Port(p = p₂′) + vol₁ = Volume(p′ = p₁′, x′ = x′, direction = -1) + vol₂ = Volume(p′ = p₂′, x′ = x′, direction = +1) + mass = Mass(f′ = (p₂′ - p₁′) * A) + flange = Flange(f = 0) end @equations begin connect(port₁, vol₁.port) @@ -132,7 +132,7 @@ end p′ end @components begin - port = Port(p=p′) + port = Port(p = p′) end @equations begin port.p ~ p′ @@ -144,20 +144,20 @@ end c = 1000 end @components begin - flange = Flange(f=0) + flange = Flange(f = 0) end @equations begin - flange.f ~ c*flange.dx + flange.f ~ c * flange.dx end end @mtkmodel System begin @components begin - res₁ = Orifice(p′=300e5) - res₂ = Orifice(p′=0) - act = Actuator(p₁′=300e5, p₂′=0) - src = Source(p′=300e5) - snk = Source(p′=0) + res₁ = Orifice(p′ = 300e5) + res₂ = Orifice(p′ = 0) + act = Actuator(p₁′ = 300e5, p₂′ = 0) + src = Source(p′ = 300e5) + snk = Source(p′ = 0) dmp = Damper() end @equations begin @@ -177,14 +177,14 @@ initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12) @test SciMLBase.successful_retcode(initsol) allinit = unknowns(sys) .=> initsol[unknowns(sys)] -prob = ODEProblem(sys, allinit, (0,0.1)) +prob = ODEProblem(sys, allinit, (0, 0.1)) sol = solve(prob, Rodas5P()) # If initialized incorrectly, then it would be InitialFailure @test sol.retcode == SciMLBase.ReturnCode.Unstable @connector Flange begin dx(t), [guess = 0] - f(t), [guess = 0, connect=Flow] + f(t), [guess = 0, connect = Flow] end @mtkmodel Mass begin @@ -202,9 +202,9 @@ end # connectors flange.dx ~ dx flange.f ~ -f - + # physics - f ~ m*D(dx) + f ~ m * D(dx) end end @@ -223,16 +223,16 @@ end # connectors flange.dx ~ dx flange.f ~ -f - + # physics - f ~ d*dx + f ~ d * dx end end @mtkmodel MassDamperSystem begin @components begin - mass = Mass(;dx=100,m=10) - damper = Damper(;d=1) + mass = Mass(; dx = 100, m = 10) + damper = Damper(; d = 1) end @equations begin connect(mass.flange, damper.flange) @@ -246,7 +246,7 @@ initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12) @test SciMLBase.successful_retcode(initsol) allinit = unknowns(sys) .=> initsol[unknowns(sys)] -prob = ODEProblem(sys, allinit, (0,0.1)) +prob = ODEProblem(sys, allinit, (0, 0.1)) sol = solve(prob, Rodas5P()) # If initialized incorrectly, then it would be InitialFailure -@test sol.retcode == SciMLBase.ReturnCode.Success \ No newline at end of file +@test sol.retcode == SciMLBase.ReturnCode.Success From 47024265d771afa554248ac5b913369f013df2ac Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 00:31:16 -0600 Subject: [PATCH 09/47] Automate tagging of the initialization system to the ODEProblem --- Project.toml | 2 +- src/systems/abstractsystem.jl | 4 ++-- src/systems/diffeqs/abstractodesystem.jl | 20 +++++++++++++++++--- test/initializationsystem.jl | 19 +++++++++++++++++++ 4 files changed, 39 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index 716b58eb03..b63b94a0f9 100644 --- a/Project.toml +++ b/Project.toml @@ -94,7 +94,7 @@ PrecompileTools = "1" RecursiveArrayTools = "2.3, 3" Reexport = "0.2, 1" RuntimeGeneratedFunctions = "0.5.9" -SciMLBase = "2.0.1" +SciMLBase = "2.27" SciMLStructures = "1.0" Serialization = "1" Setfield = "0.7, 0.8, 1" diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 3a3f2aedf6..e2f1b495f4 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -537,8 +537,8 @@ function complete(sys::AbstractSystem; split = true) if split && has_index_cache(sys) @set! sys.index_cache = IndexCache(sys) end - if isdefined(sys, :initializationsystem) - @set! sys.initializationsystem = complete(get_initializationsystem(sys); split) + if isdefined(sys, :initializesystem) && get_initializesystem(sys) !== nothing + @set! sys.initializesystem = complete(get_initializesystem(sys); split) end isdefined(sys, :complete) ? (@set! sys.complete = true) : sys end diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 42667101fe..81559b40eb 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -315,7 +315,9 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, checkbounds = false, sparsity = false, analytic = nothing, - split_idxs = nothing, + split_idxs = nothing, + initializeprob = nothing, + initializeprobmap = nothing, kwargs...) where {iip, specialize} if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEFunction`") @@ -487,6 +489,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, end @set! sys.split_idxs = split_idxs + ODEFunction{iip, specialize}(f; sys = sys, jac = _jac === nothing ? nothing : _jac, @@ -495,7 +498,9 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, jac_prototype = jac_prototype, observed = observedfun, sparsity = sparsity ? jacobian_sparsity(sys) : nothing, - analytic = analytic) + analytic = analytic, + initializeprob = initializeprob, + initializeprobmap = initializeprobmap) end """ @@ -525,6 +530,8 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) sparse = false, simplify = false, eval_module = @__MODULE__, checkbounds = false, + initializeprob = nothing, + initializeprobmap = nothing, kwargs...) where {iip} if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEFunction`") @@ -596,7 +603,9 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) sys = sys, jac = _jac === nothing ? nothing : _jac, jac_prototype = jac_prototype, - observed = observedfun) + observed = observedfun, + initializeprob = initializeprob, + initializeprobmap = initializeprobmap) end function DiffEqBase.DDEFunction(sys::AbstractODESystem, args...; kwargs...) @@ -877,10 +886,15 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; check_eqs_u0(eqs, dvs, u0; kwargs...) + initializeprob = ModelingToolkit.InitializationProblem(sys, u0map, parammap) + initializeprobmap = getu(initializeprob, unknowns(sys)) + f = constructor(sys, dvs, ps, u0; ddvs = ddvs, tgrad = tgrad, jac = jac, checkbounds = checkbounds, p = p, linenumbers = linenumbers, parallel = parallel, simplify = simplify, sparse = sparse, eval_expression = eval_expression, + initializeprob = initializeprob, + initializeprobmap = initializeprobmap, kwargs...) implicit_dae ? (f, du0, u0, p) : (f, u0, p) end diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index a92698bc1a..b0a2159867 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -181,6 +181,24 @@ prob = ODEProblem(sys, allinit, (0, 0.1)) sol = solve(prob, Rodas5P()) # If initialized incorrectly, then it would be InitialFailure @test sol.retcode == SciMLBase.ReturnCode.Unstable +SciMLBase.has_initializeprob(prob.f) + +isys = ModelingToolkit.get_initializesystem(sys) +unknowns(isys) + +initprob = ModelingToolkit.InitializationProblem(sys) +sol = solve(initprob) + +unknowns(sys) + +[sys.act.vol₁.dr] + +getter = ModelingToolkit.getu(initprob, unknowns(sys)[end-1:end]) +getter(sol) + +prob.f.initializeprobmap(initsol) + +initsol[unknowns(isys)] @connector Flange begin dx(t), [guess = 0] @@ -250,3 +268,4 @@ prob = ODEProblem(sys, allinit, (0, 0.1)) sol = solve(prob, Rodas5P()) # If initialized incorrectly, then it would be InitialFailure @test sol.retcode == SciMLBase.ReturnCode.Success + From c443ac7bae403f3c408baa51c96c022bbd81fcb7 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 02:46:48 -0600 Subject: [PATCH 10/47] Fix guess and initial condition length checking --- src/systems/abstractsystem.jl | 6 +- src/systems/diffeqs/abstractodesystem.jl | 50 +++++++++--- src/systems/nonlinear/initializesystem.jl | 6 +- src/systems/systemstructure.jl | 4 +- test/initializationsystem.jl | 96 ++++++++++++++++++----- 5 files changed, 127 insertions(+), 35 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index e2f1b495f4..a57714411e 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -2268,14 +2268,14 @@ function UnPack.unpack(sys::ModelingToolkit.AbstractSystem, ::Val{p}) where {p} end """ - missing_variable_defaults(sys::AbstractSystem, default = 0.0) + missing_variable_defaults(sys::AbstractSystem, default = 0.0; subset = unknowns(sys)) returns a `Vector{Pair}` of variables set to `default` which are missing from `get_defaults(sys)`. The `default` argument can be a single value or vector to set the missing defaults respectively. """ -function missing_variable_defaults(sys::AbstractSystem, default = 0.0) +function missing_variable_defaults(sys::AbstractSystem, default = 0.0; subset = unknowns(sys)) varmap = get_defaults(sys) varmap = Dict(Symbolics.diff2term(value(k)) => value(varmap[k]) for k in keys(varmap)) - missingvars = setdiff(unknowns(sys), keys(varmap)) + missingvars = setdiff(subset, keys(varmap)) ds = Pair[] n = length(missingvars) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 81559b40eb..1be3ab4037 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -848,18 +848,34 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; tofloat = true, symbolic_u0 = false, u0_constructor = identity, + guesses = Dict(), + warn_initialize_determined = true, kwargs...) eqs = equations(sys) dvs = unknowns(sys) ps = full_parameters(sys) iv = get_iv(sys) + initializeprob = ModelingToolkit.InitializationProblem(sys, u0map, parammap; guesses, warn_initialize_determined) + initializeprobmap = getu(initializeprob, unknowns(sys)) + + # Append zeros to the variables which are determined by the initialization system + # This essentially bypasses the check for if initial conditions are defined for DAEs + # since they will be checked in the initialization problem's construction + # TODO: make check for if a DAE cheaper than calculating the mass matrix a second time! + if implicit_dae || calculate_massmatrix(sys) !== I + zerovars = setdiff(unknowns(sys),defaults(sys)) .=> 0.0 + trueinit = identity.([zerovars;u0map]) + else + trueinit = u0map + end + if has_index_cache(sys) && get_index_cache(sys) !== nothing - u0, defs = get_u0(sys, u0map, parammap; symbolic_u0) + u0, defs = get_u0(sys, trueinit, parammap; symbolic_u0) p = MTKParameters(sys, parammap) else u0, p, defs = get_u0_p(sys, - u0map, + trueinit, parammap; tofloat, use_union, @@ -886,9 +902,6 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; check_eqs_u0(eqs, dvs, u0; kwargs...) - initializeprob = ModelingToolkit.InitializationProblem(sys, u0map, parammap) - initializeprobmap = getu(initializeprob, unknowns(sys)) - f = constructor(sys, dvs, ps, u0; ddvs = ddvs, tgrad = tgrad, jac = jac, checkbounds = checkbounds, p = p, linenumbers = linenumbers, parallel = parallel, simplify = simplify, @@ -998,13 +1011,14 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map = parammap = DiffEqBase.NullParameters(); callback = nothing, check_length = true, + warn_initialize_determined = true, kwargs...) where {iip, specialize} if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`") end f, u0, p = process_DEProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; t = tspan !== nothing ? tspan[1] : tspan, - check_length, kwargs...) + check_length, warn_initialize_determined, kwargs...) cbs = process_events(sys; callback, kwargs...) inits = [] if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing @@ -1069,13 +1083,14 @@ end function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan, parammap = DiffEqBase.NullParameters(); + warn_initialize_determined = true, check_length = true, kwargs...) where {iip} if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEProblem`") end f, du0, u0, p = process_DEProblem(DAEFunction{iip}, sys, u0map, parammap; implicit_dae = true, du0map = du0map, check_length, - kwargs...) + warn_initialize_determined, kwargs...) diffvars = collect_differential_variables(sys) sts = unknowns(sys) differential_vars = map(Base.Fix2(in, diffvars), sts) @@ -1496,18 +1511,33 @@ end function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], parammap = DiffEqBase.NullParameters(); + guesses = [], check_length = true, + warn_initialize_determined = true, kwargs...) where {iip, specialize} if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`") end - isys = get_initializesystem(sys) + if isempty(u0map) + isys = get_initializesystem(sys) + else + isys = structural_simplify(generate_initializesystem(sys; u0map); fully_determined = false) + end + neqs = length(equations(isys)) nunknown = length(unknowns(isys)) + + if warn_initialize_determined && neqs > nunknown + @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false." + end + if warn_initialize_determined && neqs < nunknown + @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false." + end + if neqs == nunknown - NonlinearProblem(isys, u0map, parammap) + NonlinearProblem(isys, guesses, parammap) else - NonlinearLeastSquaresProblem(isys, u0map, parammap) + NonlinearLeastSquaresProblem(isys, guesses, parammap) end end diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 80384bde4c..8b7ce8c9e0 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -3,7 +3,9 @@ $(TYPEDSIGNATURES) Generate `NonlinearSystem` which initializes an ODE problem from specified initial conditions of an `ODESystem`. """ -function generate_initializesystem(sys::ODESystem; name = nameof(sys), +function generate_initializesystem(sys::ODESystem; + u0map = Dict(), + name = nameof(sys), guesses = Dict(), check_defguess = false, kwargs...) sts, eqs = unknowns(sys), equations(sys) idxs_diff = isdiffeq.(eqs) @@ -13,7 +15,7 @@ function generate_initializesystem(sys::ODESystem; name = nameof(sys), # Start the equations list with algebraic equations eqs_ics = eqs[idxs_alge] u0 = Vector{Pair}(undef, 0) - defs = defaults(sys) + defs = merge(defaults(sys),todict(u0map)) full_states = [sts; getfield.((observed(sys)), :lhs)] diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index f18f5c5c89..9bba5bef5f 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -615,7 +615,7 @@ function structural_simplify!(state::TearingState, io = nothing; simplify = fals end function _structural_simplify!(state::TearingState, io; simplify = false, - check_consistency = true, fully_determined = true, warn_initialize_determined = true, + check_consistency = true, fully_determined = true, warn_initialize_determined = false, kwargs...) check_consistency &= fully_determined has_io = io !== nothing @@ -644,7 +644,7 @@ function _structural_simplify!(state::TearingState, io; simplify = false, neqs = length(equations(isys)) nunknown = length(unknowns(isys)) if warn_initialize_determined && neqs > nunknown - @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares" + @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares if $(nunknown - neqs) defaults are not supplied at construction time." end if warn_initialize_determined && neqs < nunknown @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares" diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index b0a2159867..147900310d 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -1,6 +1,49 @@ using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve, Test using ModelingToolkit: t_nounits as t, D_nounits as D +@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) + +initprob = ModelingToolkit.InitializationProblem(pend, [], [g => 1]; guesses = [ModelingToolkit.missing_variable_defaults(pend); x => 1; y => 0.2]) +conditions = getfield.(equations(initprob.f.sys),:rhs) + +@test initprob isa NonlinearLeastSquaresProblem +sol = solve(initprob) +@test SciMLBase.successful_retcode(sol) +@test maximum(abs.(sol[conditions])) < 1e-14 + +initprob = ModelingToolkit.InitializationProblem(pend, [x => 1, y => 0], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) +@test initprob isa NonlinearProblem +sol = solve(initprob) +@test SciMLBase.successful_retcode(sol) +@test sol.u == [1.0,0.0,0.0,0.0] +@test maximum(abs.(sol[conditions])) < 1e-14 + +initprob = ModelingToolkit.InitializationProblem(pend, [], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) +@test initprob isa NonlinearLeastSquaresProblem +sol = solve(initprob) +@test !SciMLBase.successful_retcode(sol) + +prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], guesses = ModelingToolkit.missing_variable_defaults(pend)) +prob.f.initializeprob isa NonlinearProblem +sol = solve(prob.f.initializeprob) +@test maximum(abs.(sol[conditions])) < 1e-14 +sol = solve(prob, Rodas5P()) +@test maximum(abs.(sol[conditions][1])) < 1e-14 + +prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = ModelingToolkit.missing_variable_defaults(pend)) +prob.f.initializeprob isa NonlinearLeastSquaresProblem +sol = solve(prob.f.initializeprob) +@test maximum(abs.(sol[conditions])) < 1e-14 +sol = solve(prob, Rodas5P()) +@test maximum(abs.(sol[conditions][1])) < 1e-14 + @connector Port begin p(t) dm(t) = 0, [connect = Flow] @@ -171,34 +214,26 @@ end @mtkbuild sys = System() initprob = ModelingToolkit.InitializationProblem(sys) +conditions = getfield.(equations(initprob.f.sys),:rhs) + @test initprob isa NonlinearLeastSquaresProblem @test length(initprob.u0) == 2 initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12) @test SciMLBase.successful_retcode(initsol) +@test maximum(abs.(initsol[conditions])) < 1e-14 allinit = unknowns(sys) .=> initsol[unknowns(sys)] prob = ODEProblem(sys, allinit, (0, 0.1)) -sol = solve(prob, Rodas5P()) +sol = solve(prob, Rodas5P(), initializealg = BrownFullBasicInit()) # If initialized incorrectly, then it would be InitialFailure @test sol.retcode == SciMLBase.ReturnCode.Unstable -SciMLBase.has_initializeprob(prob.f) - -isys = ModelingToolkit.get_initializesystem(sys) -unknowns(isys) - -initprob = ModelingToolkit.InitializationProblem(sys) -sol = solve(initprob) - -unknowns(sys) - -[sys.act.vol₁.dr] - -getter = ModelingToolkit.getu(initprob, unknowns(sys)[end-1:end]) -getter(sol) +@test maximum(abs.(initsol[conditions][1])) < 1e-14 -prob.f.initializeprobmap(initsol) - -initsol[unknowns(isys)] +prob = ODEProblem(sys, [], (0, 0.1), check=false) +sol = solve(prob, Rodas5P()) +# If initialized incorrectly, then it would be InitialFailure +@test sol.retcode == SciMLBase.ReturnCode.Unstable +@test maximum(abs.(initsol[conditions][1])) < 1e-14 @connector Flange begin dx(t), [guess = 0] @@ -269,3 +304,28 @@ sol = solve(prob, Rodas5P()) # If initialized incorrectly, then it would be InitialFailure @test sol.retcode == SciMLBase.ReturnCode.Success +prob = ODEProblem(sys, [], (0, 0.1)) +sol = solve(prob, Rodas5P()) +@test sol.retcode == SciMLBase.ReturnCode.Success + +### Ensure that non-DAEs still throw for missing variables without the initialize system + +@parameters σ ρ β +@variables x(t) y(t) z(t) + +eqs = [D(D(x)) ~ σ * (y - x), + D(y) ~ x * (ρ - z) - y, + D(z) ~ x * y - β * z] + +@mtkbuild sys = ODESystem(eqs, t) + +u0 = [D(x) => 2.0, + y => 0.0, + z => 0.0] + +p = [σ => 28.0, + ρ => 10.0, + β => 8 / 3] + +tspan = (0.0, 100.0) +@test_throws ArgumentError prob = ODEProblem(sys, u0, tspan, p, jac = true) \ No newline at end of file From 67f6f24b9288ab838911341d2244783f29aa2cae Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 09:34:05 -0600 Subject: [PATCH 11/47] Bump ordinarydiffeq --- Project.toml | 4 ++-- src/systems/diffeqs/abstractodesystem.jl | 7 ++++--- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index b63b94a0f9..76ccd0e386 100644 --- a/Project.toml +++ b/Project.toml @@ -89,12 +89,12 @@ Libdl = "1" LinearAlgebra = "1" MLStyle = "0.4.17" NaNMath = "0.3, 1" -OrdinaryDiffEq = "6" +OrdinaryDiffEq = "6.72.0" PrecompileTools = "1" RecursiveArrayTools = "2.3, 3" Reexport = "0.2, 1" RuntimeGeneratedFunctions = "0.5.9" -SciMLBase = "2.27" +SciMLBase = "2.28.0" SciMLStructures = "1.0" Serialization = "1" Setfield = "0.7, 0.8, 1" diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 1be3ab4037..fb1b5d3fb4 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -856,17 +856,18 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; ps = full_parameters(sys) iv = get_iv(sys) - initializeprob = ModelingToolkit.InitializationProblem(sys, u0map, parammap; guesses, warn_initialize_determined) - initializeprobmap = getu(initializeprob, unknowns(sys)) - # Append zeros to the variables which are determined by the initialization system # This essentially bypasses the check for if initial conditions are defined for DAEs # since they will be checked in the initialization problem's construction # TODO: make check for if a DAE cheaper than calculating the mass matrix a second time! if implicit_dae || calculate_massmatrix(sys) !== I + initializeprob = ModelingToolkit.InitializationProblem(sys, u0map, parammap; guesses, warn_initialize_determined) + initializeprobmap = getu(initializeprob, unknowns(sys)) zerovars = setdiff(unknowns(sys),defaults(sys)) .=> 0.0 trueinit = identity.([zerovars;u0map]) else + initializeprob = nothing + initializeprobmap = nothing trueinit = u0map end From 7c0c423b9c7204f92329ed8ba08de6206991d2d3 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 10:19:44 -0600 Subject: [PATCH 12/47] format --- src/systems/abstractsystem.jl | 3 ++- src/systems/diffeqs/abstractodesystem.jl | 20 +++++++------- src/systems/nonlinear/initializesystem.jl | 4 +-- test/initializationsystem.jl | 33 ++++++++++++----------- 4 files changed, 33 insertions(+), 27 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index a57714411e..3ff6527755 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -2272,7 +2272,8 @@ end returns a `Vector{Pair}` of variables set to `default` which are missing from `get_defaults(sys)`. The `default` argument can be a single value or vector to set the missing defaults respectively. """ -function missing_variable_defaults(sys::AbstractSystem, default = 0.0; subset = unknowns(sys)) +function missing_variable_defaults( + sys::AbstractSystem, default = 0.0; subset = unknowns(sys)) varmap = get_defaults(sys) varmap = Dict(Symbolics.diff2term(value(k)) => value(varmap[k]) for k in keys(varmap)) missingvars = setdiff(subset, keys(varmap)) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index fb1b5d3fb4..fd99a0510b 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -315,8 +315,8 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, checkbounds = false, sparsity = false, analytic = nothing, - split_idxs = nothing, - initializeprob = nothing, + split_idxs = nothing, + initializeprob = nothing, initializeprobmap = nothing, kwargs...) where {iip, specialize} if !iscomplete(sys) @@ -530,7 +530,7 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) sparse = false, simplify = false, eval_module = @__MODULE__, checkbounds = false, - initializeprob = nothing, + initializeprob = nothing, initializeprobmap = nothing, kwargs...) where {iip} if !iscomplete(sys) @@ -861,12 +861,13 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; # since they will be checked in the initialization problem's construction # TODO: make check for if a DAE cheaper than calculating the mass matrix a second time! if implicit_dae || calculate_massmatrix(sys) !== I - initializeprob = ModelingToolkit.InitializationProblem(sys, u0map, parammap; guesses, warn_initialize_determined) + initializeprob = ModelingToolkit.InitializationProblem( + sys, u0map, parammap; guesses, warn_initialize_determined) initializeprobmap = getu(initializeprob, unknowns(sys)) - zerovars = setdiff(unknowns(sys),defaults(sys)) .=> 0.0 - trueinit = identity.([zerovars;u0map]) + zerovars = setdiff(unknowns(sys), defaults(sys)) .=> 0.0 + trueinit = identity.([zerovars; u0map]) else - initializeprob = nothing + initializeprob = nothing initializeprobmap = nothing trueinit = u0map end @@ -907,7 +908,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; checkbounds = checkbounds, p = p, linenumbers = linenumbers, parallel = parallel, simplify = simplify, sparse = sparse, eval_expression = eval_expression, - initializeprob = initializeprob, + initializeprob = initializeprob, initializeprobmap = initializeprobmap, kwargs...) implicit_dae ? (f, du0, u0, p) : (f, u0, p) @@ -1523,7 +1524,8 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0map = if isempty(u0map) isys = get_initializesystem(sys) else - isys = structural_simplify(generate_initializesystem(sys; u0map); fully_determined = false) + isys = structural_simplify( + generate_initializesystem(sys; u0map); fully_determined = false) end neqs = length(equations(isys)) diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 8b7ce8c9e0..51e9c49480 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -3,7 +3,7 @@ $(TYPEDSIGNATURES) Generate `NonlinearSystem` which initializes an ODE problem from specified initial conditions of an `ODESystem`. """ -function generate_initializesystem(sys::ODESystem; +function generate_initializesystem(sys::ODESystem; u0map = Dict(), name = nameof(sys), guesses = Dict(), check_defguess = false, kwargs...) @@ -15,7 +15,7 @@ function generate_initializesystem(sys::ODESystem; # Start the equations list with algebraic equations eqs_ics = eqs[idxs_alge] u0 = Vector{Pair}(undef, 0) - defs = merge(defaults(sys),todict(u0map)) + defs = merge(defaults(sys), todict(u0map)) full_states = [sts; getfield.((observed(sys)), :lhs)] diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index 147900310d..2e6fc8c95c 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -3,41 +3,44 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @parameters g @variables x(t) y(t) [state_priority = 10] λ(t) -eqs = [ - D(D(x)) ~ λ * x +eqs = [D(D(x)) ~ λ * x D(D(y)) ~ λ * y - g - x^2 + y^2 ~ 1 - ] -@mtkbuild pend = ODESystem(eqs,t) + x^2 + y^2 ~ 1] +@mtkbuild pend = ODESystem(eqs, t) -initprob = ModelingToolkit.InitializationProblem(pend, [], [g => 1]; guesses = [ModelingToolkit.missing_variable_defaults(pend); x => 1; y => 0.2]) -conditions = getfield.(equations(initprob.f.sys),:rhs) +initprob = ModelingToolkit.InitializationProblem(pend, [], [g => 1]; + guesses = [ModelingToolkit.missing_variable_defaults(pend); x => 1; y => 0.2]) +conditions = getfield.(equations(initprob.f.sys), :rhs) @test initprob isa NonlinearLeastSquaresProblem sol = solve(initprob) @test SciMLBase.successful_retcode(sol) @test maximum(abs.(sol[conditions])) < 1e-14 -initprob = ModelingToolkit.InitializationProblem(pend, [x => 1, y => 0], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) +initprob = ModelingToolkit.InitializationProblem(pend, [x => 1, y => 0], [g => 1]; + guesses = ModelingToolkit.missing_variable_defaults(pend)) @test initprob isa NonlinearProblem sol = solve(initprob) @test SciMLBase.successful_retcode(sol) -@test sol.u == [1.0,0.0,0.0,0.0] +@test sol.u == [1.0, 0.0, 0.0, 0.0] @test maximum(abs.(sol[conditions])) < 1e-14 -initprob = ModelingToolkit.InitializationProblem(pend, [], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) +initprob = ModelingToolkit.InitializationProblem( + pend, [], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) @test initprob isa NonlinearLeastSquaresProblem sol = solve(initprob) @test !SciMLBase.successful_retcode(sol) -prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], guesses = ModelingToolkit.missing_variable_defaults(pend)) +prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], + guesses = ModelingToolkit.missing_variable_defaults(pend)) prob.f.initializeprob isa NonlinearProblem sol = solve(prob.f.initializeprob) @test maximum(abs.(sol[conditions])) < 1e-14 sol = solve(prob, Rodas5P()) @test maximum(abs.(sol[conditions][1])) < 1e-14 -prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = ModelingToolkit.missing_variable_defaults(pend)) +prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], + guesses = ModelingToolkit.missing_variable_defaults(pend)) prob.f.initializeprob isa NonlinearLeastSquaresProblem sol = solve(prob.f.initializeprob) @test maximum(abs.(sol[conditions])) < 1e-14 @@ -214,7 +217,7 @@ end @mtkbuild sys = System() initprob = ModelingToolkit.InitializationProblem(sys) -conditions = getfield.(equations(initprob.f.sys),:rhs) +conditions = getfield.(equations(initprob.f.sys), :rhs) @test initprob isa NonlinearLeastSquaresProblem @test length(initprob.u0) == 2 @@ -229,7 +232,7 @@ sol = solve(prob, Rodas5P(), initializealg = BrownFullBasicInit()) @test sol.retcode == SciMLBase.ReturnCode.Unstable @test maximum(abs.(initsol[conditions][1])) < 1e-14 -prob = ODEProblem(sys, [], (0, 0.1), check=false) +prob = ODEProblem(sys, [], (0, 0.1), check = false) sol = solve(prob, Rodas5P()) # If initialized incorrectly, then it would be InitialFailure @test sol.retcode == SciMLBase.ReturnCode.Unstable @@ -328,4 +331,4 @@ p = [σ => 28.0, β => 8 / 3] tspan = (0.0, 100.0) -@test_throws ArgumentError prob = ODEProblem(sys, u0, tspan, p, jac = true) \ No newline at end of file +@test_throws ArgumentError prob=ODEProblem(sys, u0, tspan, p, jac = true) From fb4b987df34a8e6eaed57582d1bd67d0b3aec366 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 11:58:36 -0600 Subject: [PATCH 13/47] stop initialize on clocks --- src/systems/diffeqs/abstractodesystem.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index fd99a0510b..584bbcc62a 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -860,7 +860,9 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; # This essentially bypasses the check for if initial conditions are defined for DAEs # since they will be checked in the initialization problem's construction # TODO: make check for if a DAE cheaper than calculating the mass matrix a second time! - if implicit_dae || calculate_massmatrix(sys) !== I + ci = infer_clocks!(ClockInference(TearingState(sys))) + # TODO: make it work with clocks + if (implicit_dae || calculate_massmatrix(sys) !== I) && all(isequal(Continuous()),ci.var_domain) initializeprob = ModelingToolkit.InitializationProblem( sys, u0map, parammap; guesses, warn_initialize_determined) initializeprobmap = getu(initializeprob, unknowns(sys)) From 16ce722cd12fb167fb716ea559232b3d309a7a81 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 13:01:47 -0600 Subject: [PATCH 14/47] format --- src/systems/diffeqs/abstractodesystem.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 584bbcc62a..d21f4e7766 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -862,7 +862,8 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; # TODO: make check for if a DAE cheaper than calculating the mass matrix a second time! ci = infer_clocks!(ClockInference(TearingState(sys))) # TODO: make it work with clocks - if (implicit_dae || calculate_massmatrix(sys) !== I) && all(isequal(Continuous()),ci.var_domain) + if (implicit_dae || calculate_massmatrix(sys) !== I) && + all(isequal(Continuous()), ci.var_domain) initializeprob = ModelingToolkit.InitializationProblem( sys, u0map, parammap; guesses, warn_initialize_determined) initializeprobmap = getu(initializeprob, unknowns(sys)) From 7a844440557c77e5135e3b2161905933def6b8ee Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 13:53:26 -0600 Subject: [PATCH 15/47] fix direct structural simplification for initialization systems --- src/systems/systemstructure.jl | 4 +++- test/clock.jl | 1 + 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 9bba5bef5f..aa0348a1d9 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -636,7 +636,9 @@ function _structural_simplify!(state::TearingState, io; simplify = false, fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)] @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns) - if sys isa ODESystem + ci = infer_clocks!(ClockInference(state)) + # TODO: make it work with clocks + if sys isa ODESystem && all(isequal(Continuous()), ci.var_domain) && !all(all(x->!(typeof(x) <: Union{Sample,Hold,ShiftIndex}),io)) isys = ModelingToolkit.generate_initializesystem(sys) !isempty(equations(isys)) && (isys = structural_simplify(isys; fully_determined = false)) diff --git a/test/clock.jl b/test/clock.jl index 26b416fe26..0ef3718f5d 100644 --- a/test/clock.jl +++ b/test/clock.jl @@ -62,6 +62,7 @@ By inference: => Shift(x, 0, dt) := (Shift(x, -1, dt) + dt) / (1 - dt) # Discrete system =# + ci, varmap = infer_clocks(sys) eqmap = ci.eq_domain tss, inputs = ModelingToolkit.split_system(deepcopy(ci)) From ae4be0f68372b5c4d9ae7ba3533aa3eb71e97f27 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 13:57:18 -0600 Subject: [PATCH 16/47] format --- src/systems/systemstructure.jl | 3 ++- test/clock.jl | 1 - 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index aa0348a1d9..f4d5bad2d3 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -638,7 +638,8 @@ function _structural_simplify!(state::TearingState, io; simplify = false, ci = infer_clocks!(ClockInference(state)) # TODO: make it work with clocks - if sys isa ODESystem && all(isequal(Continuous()), ci.var_domain) && !all(all(x->!(typeof(x) <: Union{Sample,Hold,ShiftIndex}),io)) + if sys isa ODESystem && all(isequal(Continuous()), ci.var_domain) && + !all(all(x -> !(typeof(x) <: Union{Sample, Hold, ShiftIndex}), io)) isys = ModelingToolkit.generate_initializesystem(sys) !isempty(equations(isys)) && (isys = structural_simplify(isys; fully_determined = false)) diff --git a/test/clock.jl b/test/clock.jl index 0ef3718f5d..26b416fe26 100644 --- a/test/clock.jl +++ b/test/clock.jl @@ -62,7 +62,6 @@ By inference: => Shift(x, 0, dt) := (Shift(x, -1, dt) + dt) / (1 - dt) # Discrete system =# - ci, varmap = infer_clocks(sys) eqmap = ci.eq_domain tss, inputs = ModelingToolkit.split_system(deepcopy(ci)) From dcffe8f7ec1292075bb5eb0c1f1a0ba15fb4a421 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 14:19:53 -0600 Subject: [PATCH 17/47] check for io first --- src/systems/systemstructure.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index f4d5bad2d3..bee7d97a51 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -638,8 +638,8 @@ function _structural_simplify!(state::TearingState, io; simplify = false, ci = infer_clocks!(ClockInference(state)) # TODO: make it work with clocks - if sys isa ODESystem && all(isequal(Continuous()), ci.var_domain) && - !all(all(x -> !(typeof(x) <: Union{Sample, Hold, ShiftIndex}), io)) + if sys isa ODESystem && all(isequal(Continuous()), ci.var_domain) && (!has_io || + !all(all(x -> !(typeof(x) <: Union{Sample, Hold, ShiftIndex}), io))) isys = ModelingToolkit.generate_initializesystem(sys) !isempty(equations(isys)) && (isys = structural_simplify(isys; fully_determined = false)) From 78b58329eae8229f5bd7cd2dc39d5f71a3153a55 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 14:23:35 -0600 Subject: [PATCH 18/47] format --- src/systems/systemstructure.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index bee7d97a51..3f022640fe 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -638,8 +638,9 @@ function _structural_simplify!(state::TearingState, io; simplify = false, ci = infer_clocks!(ClockInference(state)) # TODO: make it work with clocks - if sys isa ODESystem && all(isequal(Continuous()), ci.var_domain) && (!has_io || - !all(all(x -> !(typeof(x) <: Union{Sample, Hold, ShiftIndex}), io))) + if sys isa ODESystem && all(isequal(Continuous()), ci.var_domain) && + (!has_io || + !all(all(x -> !(typeof(x) <: Union{Sample, Hold, ShiftIndex}), io))) isys = ModelingToolkit.generate_initializesystem(sys) !isempty(equations(isys)) && (isys = structural_simplify(isys; fully_determined = false)) From 6e336ba0ac2700dc60481438146e5ba373250116 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 17:50:49 -0600 Subject: [PATCH 19/47] Properly drop from defaults --- src/systems/diffeqs/abstractodesystem.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index d21f4e7766..7303c74596 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -867,7 +867,8 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; initializeprob = ModelingToolkit.InitializationProblem( sys, u0map, parammap; guesses, warn_initialize_determined) initializeprobmap = getu(initializeprob, unknowns(sys)) - zerovars = setdiff(unknowns(sys), defaults(sys)) .=> 0.0 + + zerovars = setdiff(unknowns(sys), keys(defaults(sys))) .=> 0.0 trueinit = identity.([zerovars; u0map]) else initializeprob = nothing From dfee0efea0342e37d660979384402a3f83f53130 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 19:27:17 -0600 Subject: [PATCH 20/47] handle the Vector{Float} case --- src/systems/diffeqs/abstractodesystem.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 7303c74596..92c53ecb11 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -864,6 +864,9 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; # TODO: make it work with clocks if (implicit_dae || calculate_massmatrix(sys) !== I) && all(isequal(Continuous()), ci.var_domain) + if eltype(u0map) <: Number + u0map = unknowns(sys) .=> u0map + end initializeprob = ModelingToolkit.InitializationProblem( sys, u0map, parammap; guesses, warn_initialize_determined) initializeprobmap = getu(initializeprob, unknowns(sys)) From f62163a1ffefbd645935229c7c0abf84574339d8 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 19:57:01 -0600 Subject: [PATCH 21/47] handle nothing equations case --- src/systems/diffeqs/abstractodesystem.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 92c53ecb11..6fbca2dc4b 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1535,6 +1535,9 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0map = generate_initializesystem(sys; u0map); fully_determined = false) end + if equations(isys) === nothing + return NonlinearProblem(isys, guesses, parammap) + end neqs = length(equations(isys)) nunknown = length(unknowns(isys)) From 80ece19da6e98edd68e5db90492de28db92a9033 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 20:37:18 -0600 Subject: [PATCH 22/47] handle no structural simplify --- src/systems/diffeqs/abstractodesystem.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 6fbca2dc4b..02f5ed868b 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1528,16 +1528,15 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0map = error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`") end - if isempty(u0map) + if isempty(u0map) && get_initializesystem(sys) !== nothing isys = get_initializesystem(sys) + elseif isempty(u0map) && get_initializesystem(sys) === nothing + isys = structural_simplify(generate_initializesystem(sys); fully_determined = false) else isys = structural_simplify( generate_initializesystem(sys; u0map); fully_determined = false) end - if equations(isys) === nothing - return NonlinearProblem(isys, guesses, parammap) - end neqs = length(equations(isys)) nunknown = length(unknowns(isys)) From 97c1d9df1571d144b5c97c1c8691d991ce328807 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 26 Feb 2024 20:59:26 -0600 Subject: [PATCH 23/47] new error --- test/odesystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/odesystem.jl b/test/odesystem.jl index efea02b628..fa24ab20ac 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -300,7 +300,7 @@ prob3 = ODEProblem(sys, u0, tspan, p, jac = true, sparse = true) #SparseMatrixCS @test prob3.f.jac_prototype isa SparseMatrixCSC prob3 = ODEProblem(sys, u0, tspan, p, jac = true, sparsity = true) @test prob3.f.sparsity isa SparseMatrixCSC -@test_throws ArgumentError ODEProblem(sys, zeros(5), tspan, p) +@test_throws DimensionMismatch ODEProblem(sys, zeros(5), tspan, p) for (prob, atol) in [(prob1, 1e-12), (prob2, 1e-12), (prob3, 1e-12)] local sol sol = solve(prob, Rodas5()) From 8521ddb427cdd1918fb04af9aa8290e6822f3c66 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 01:06:33 -0600 Subject: [PATCH 24/47] Remove scalarization in NonlinearSystem --- src/systems/nonlinear/nonlinearsystem.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index 798ed10e8d..9459ba0a9b 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -142,7 +142,6 @@ function NonlinearSystem(eqs, unknowns, ps; defaults = todict(defaults) defaults = Dict{Any, Any}(value(k) => value(v) for (k, v) in pairs(defaults)) - unknowns = scalarize(unknowns) unknowns, ps = value.(unknowns), value.(ps) var_to_name = Dict() process_variables!(var_to_name, defaults, unknowns) @@ -362,9 +361,11 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem, u0map, para dvs = unknowns(sys) ps = parameters(sys) - u0, p, defs = get_u0_p(sys, u0map, parammap; tofloat, use_union) if has_index_cache(sys) && get_index_cache(sys) !== nothing + u0, defs = get_u0(sys, u0map, parammap) p = MTKParameters(sys, parammap) + else + u0, p, defs = get_u0_p(sys, u0map, parammap; tofloat, use_union) end check_eqs_u0(eqs, dvs, u0; kwargs...) From 2efc30aeb94c4121ab0df148908dae6ff8851abd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 01:56:27 -0600 Subject: [PATCH 25/47] Try other tests --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 605825bc0b..efdbda9607 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,7 +30,7 @@ end @safetestset "Dynamic Quantities Test" include("dq_units.jl") @safetestset "Unitful Quantities Test" include("units.jl") @safetestset "LabelledArrays Test" include("labelledarrays.jl") - @safetestset "Mass Matrix Test" include("mass_matrix.jl") + #@safetestset "Mass Matrix Test" include("mass_matrix.jl") @safetestset "SteadyStateSystem Test" include("steadystatesystems.jl") @safetestset "SDESystem Test" include("sdesystem.jl") @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl") From 9df0b4b4585e70d25a9c5f1913d6888de6d12ee0 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 03:31:53 -0600 Subject: [PATCH 26/47] remove old initialize --- test/nonlinearsystem.jl | 43 ----------------------------------------- 1 file changed, 43 deletions(-) diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl index f8e4541a75..49de424fa4 100644 --- a/test/nonlinearsystem.jl +++ b/test/nonlinearsystem.jl @@ -242,46 +242,3 @@ testdict = Dict([:test => 1]) @test prob_.u0 == [1.0, 2.0, 1.0] @test prob_.p == MTKParameters(sys, [a => 2.0, b => 1.0, c => 1.0]) end - -@testset "Initialization System" begin - # Define the Lotka Volterra system which begins at steady state - @parameters t - pars = @parameters a=1.5 b=1.0 c=3.0 d=1.0 dx_ss=1e-5 - - vars = @variables begin - dx(t), - dy(t), - (x(t) = dx ~ dx_ss), [guess = 0.5] - (y(t) = dy ~ 0), [guess = -0.5] - end - - D = Differential(t) - - eqs = [dx ~ a * x - b * x * y - dy ~ -c * y + d * x * y - D(x) ~ dx - D(y) ~ dy] - - @named sys = ODESystem(eqs, t, vars, pars) - - sys_simple = structural_simplify(sys) - - # Set up the initialization system - sys_init = initializesystem(sys_simple) - - sys_init_simple = structural_simplify(sys_init) - - prob = NonlinearProblem(sys_init_simple, - get_default_or_guess.(unknowns(sys_init_simple))) - - @test prob.u0 == [0.5, -0.5] - - sol = solve(prob) - @test sol.retcode == SciMLBase.ReturnCode.Success - - # Confirm for all the unknowns of the non-simplified system - @test all(.≈(sol[unknowns(sys)], [1e-5, 0, 1e-5 / 1.5, 0]; atol = 1e-8)) - - # Confirm for all the unknowns of the simplified system - @test all(.≈(sol[unknowns(sys_simple)], [1e-5 / 1.5, 0]; atol = 1e-8)) -end From 3d41e04c4044cc98e52a497727badda28ce04027 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 05:21:34 -0600 Subject: [PATCH 27/47] see if the non array tests all pass --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index efdbda9607..b219edf30c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -38,7 +38,7 @@ end @safetestset "PDE Construction Test" include("pde.jl") @safetestset "JumpSystem Test" include("jumpsystem.jl") @safetestset "Constraints Test" include("constraints.jl") - @safetestset "Reduction Test" include("reduction.jl") + #@safetestset "Reduction Test" include("reduction.jl") @safetestset "Split Parameters Test" include("split_parameters.jl") @safetestset "StaticArrays Test" include("static_arrays.jl") @safetestset "Components Test" include("components.jl") From 07638741f425ed17a412245eb32916d6b40dcb3b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 05:23:43 -0600 Subject: [PATCH 28/47] fix reduction test --- test/reduction.jl | 2 +- test/runtests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/reduction.jl b/test/reduction.jl index 7a8979ba23..0dbd032447 100644 --- a/test/reduction.jl +++ b/test/reduction.jl @@ -176,7 +176,7 @@ N = 5 @variables xs[1:N] A = reshape(1:(N^2), N, N) eqs = xs .~ A * xs -@named sys′ = NonlinearSystem(eqs, xs, []) +@named sys′ = NonlinearSystem(collect(eqs), [xs], []) sys = structural_simplify(sys′) # issue 958 diff --git a/test/runtests.jl b/test/runtests.jl index b219edf30c..efdbda9607 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -38,7 +38,7 @@ end @safetestset "PDE Construction Test" include("pde.jl") @safetestset "JumpSystem Test" include("jumpsystem.jl") @safetestset "Constraints Test" include("constraints.jl") - #@safetestset "Reduction Test" include("reduction.jl") + @safetestset "Reduction Test" include("reduction.jl") @safetestset "Split Parameters Test" include("split_parameters.jl") @safetestset "StaticArrays Test" include("static_arrays.jl") @safetestset "Components Test" include("components.jl") From deaa8247520854bcc35eff057296c07d2aac146f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 05:26:15 -0600 Subject: [PATCH 29/47] better fix --- test/reduction.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/reduction.jl b/test/reduction.jl index 0dbd032447..5c86ed0ef9 100644 --- a/test/reduction.jl +++ b/test/reduction.jl @@ -175,8 +175,8 @@ nlprob.f(residual, reducedsol.u, pp) N = 5 @variables xs[1:N] A = reshape(1:(N^2), N, N) -eqs = xs .~ A * xs -@named sys′ = NonlinearSystem(collect(eqs), [xs], []) +eqs = xs ~ A * xs +@named sys′ = NonlinearSystem(eqs, [xs], []) sys = structural_simplify(sys′) # issue 958 From 456478f1867edaf76e9c9d5283fce08437491d59 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 05:41:57 -0600 Subject: [PATCH 30/47] only build initialization if simplified --- src/systems/diffeqs/abstractodesystem.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 02f5ed868b..584f1c6ff1 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -862,8 +862,11 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; # TODO: make check for if a DAE cheaper than calculating the mass matrix a second time! ci = infer_clocks!(ClockInference(TearingState(sys))) # TODO: make it work with clocks + # ModelingToolkit.get_tearing_state(sys) !== nothing => Requires structural_simplify first if (implicit_dae || calculate_massmatrix(sys) !== I) && - all(isequal(Continuous()), ci.var_domain) + all(isequal(Continuous()), ci.var_domain) && + ModelingToolkit.get_tearing_state(sys) !== nothing + if eltype(u0map) <: Number u0map = unknowns(sys) .=> u0map end From a5e275d23cdc86ab636bde676da44b2cfd61c19b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 05:43:20 -0600 Subject: [PATCH 31/47] reenable mass matrix test --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index efdbda9607..605825bc0b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,7 +30,7 @@ end @safetestset "Dynamic Quantities Test" include("dq_units.jl") @safetestset "Unitful Quantities Test" include("units.jl") @safetestset "LabelledArrays Test" include("labelledarrays.jl") - #@safetestset "Mass Matrix Test" include("mass_matrix.jl") + @safetestset "Mass Matrix Test" include("mass_matrix.jl") @safetestset "SteadyStateSystem Test" include("steadystatesystems.jl") @safetestset "SDESystem Test" include("sdesystem.jl") @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl") From f85b21905cc0976f9d65ef98550b2ad89cff0305 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 07:00:59 -0600 Subject: [PATCH 32/47] format --- src/systems/diffeqs/abstractodesystem.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 584f1c6ff1..430e9e5a85 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -866,7 +866,6 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; if (implicit_dae || calculate_massmatrix(sys) !== I) && all(isequal(Continuous()), ci.var_domain) && ModelingToolkit.get_tearing_state(sys) !== nothing - if eltype(u0map) <: Number u0map = unknowns(sys) .=> u0map end From f2559ab32fc448e07116e9ea3a77592d96b3b6db Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 07:26:23 -0600 Subject: [PATCH 33/47] Update test/odesystem.jl --- test/odesystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/odesystem.jl b/test/odesystem.jl index fa24ab20ac..efea02b628 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -300,7 +300,7 @@ prob3 = ODEProblem(sys, u0, tspan, p, jac = true, sparse = true) #SparseMatrixCS @test prob3.f.jac_prototype isa SparseMatrixCSC prob3 = ODEProblem(sys, u0, tspan, p, jac = true, sparsity = true) @test prob3.f.sparsity isa SparseMatrixCSC -@test_throws DimensionMismatch ODEProblem(sys, zeros(5), tspan, p) +@test_throws ArgumentError ODEProblem(sys, zeros(5), tspan, p) for (prob, atol) in [(prob1, 1e-12), (prob2, 1e-12), (prob3, 1e-12)] local sol sol = solve(prob, Rodas5()) From 7970cc1a672727db475cacc21323c5f1ac2b806f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 08:03:18 -0600 Subject: [PATCH 34/47] let all tests run --- test/runtests.jl | 100 ++++++++++++++++++++++++----------------------- 1 file changed, 51 insertions(+), 49 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 605825bc0b..371834f351 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,55 +16,57 @@ end @time begin if GROUP == "All" || GROUP == "InterfaceI" - @safetestset "Linear Algebra Test" include("linalg.jl") - @safetestset "AbstractSystem Test" include("abstractsystem.jl") - @safetestset "Variable Scope Tests" include("variable_scope.jl") - @safetestset "Symbolic Parameters Test" include("symbolic_parameters.jl") - @safetestset "Parsing Test" include("variable_parsing.jl") - @safetestset "Simplify Test" include("simplify.jl") - @safetestset "Direct Usage Test" include("direct.jl") - @safetestset "System Linearity Test" include("linearity.jl") - @safetestset "Input Output Test" include("input_output_handling.jl") - @safetestset "Clock Test" include("clock.jl") - @safetestset "ODESystem Test" include("odesystem.jl") - @safetestset "Dynamic Quantities Test" include("dq_units.jl") - @safetestset "Unitful Quantities Test" include("units.jl") - @safetestset "LabelledArrays Test" include("labelledarrays.jl") - @safetestset "Mass Matrix Test" include("mass_matrix.jl") - @safetestset "SteadyStateSystem Test" include("steadystatesystems.jl") - @safetestset "SDESystem Test" include("sdesystem.jl") - @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl") - @safetestset "InitializationSystem Test" include("initializationsystem.jl") - @safetestset "PDE Construction Test" include("pde.jl") - @safetestset "JumpSystem Test" include("jumpsystem.jl") - @safetestset "Constraints Test" include("constraints.jl") - @safetestset "Reduction Test" include("reduction.jl") - @safetestset "Split Parameters Test" include("split_parameters.jl") - @safetestset "StaticArrays Test" include("static_arrays.jl") - @safetestset "Components Test" include("components.jl") - @safetestset "Model Parsing Test" include("model_parsing.jl") - @safetestset "print_tree" include("print_tree.jl") - @safetestset "Error Handling" include("error_handling.jl") - @safetestset "StructuralTransformations" include("structural_transformation/runtests.jl") - @safetestset "State Selection Test" include("state_selection.jl") - @safetestset "Symbolic Event Test" include("symbolic_events.jl") - @safetestset "Stream Connect Test" include("stream_connectors.jl") - @safetestset "Domain Connect Test" include("domain_connectors.jl") - @safetestset "Lowering Integration Test" include("lowering_solving.jl") - @safetestset "Test Big System Usage" include("bigsystem.jl") - @safetestset "Dependency Graph Test" include("dep_graphs.jl") - @safetestset "Function Registration Test" include("function_registration.jl") - @safetestset "Precompiled Modules Test" include("precompile_test.jl") - @safetestset "Variable Utils Test" include("variable_utils.jl") - @safetestset "Variable Metadata Test" include("test_variable_metadata.jl") - @safetestset "DAE Jacobians Test" include("dae_jacobian.jl") - @safetestset "Jacobian Sparsity" include("jacobiansparsity.jl") - @safetestset "Modelingtoolkitize Test" include("modelingtoolkitize.jl") - @safetestset "OptimizationSystem Test" include("optimizationsystem.jl") - @safetestset "FuncAffect Test" include("funcaffect.jl") - @safetestset "Constants Test" include("constants.jl") - @safetestset "Parameter Dependency Test" include("parameter_dependencies.jl") - @safetestset "Generate Custom Function Test" include("generate_custom_function.jl") + @testset "InterfaceI" begin + @safetestset "Linear Algebra Test" include("linalg.jl") + @safetestset "AbstractSystem Test" include("abstractsystem.jl") + @safetestset "Variable Scope Tests" include("variable_scope.jl") + @safetestset "Symbolic Parameters Test" include("symbolic_parameters.jl") + @safetestset "Parsing Test" include("variable_parsing.jl") + @safetestset "Simplify Test" include("simplify.jl") + @safetestset "Direct Usage Test" include("direct.jl") + @safetestset "System Linearity Test" include("linearity.jl") + @safetestset "Input Output Test" include("input_output_handling.jl") + @safetestset "Clock Test" include("clock.jl") + @safetestset "ODESystem Test" include("odesystem.jl") + @safetestset "Dynamic Quantities Test" include("dq_units.jl") + @safetestset "Unitful Quantities Test" include("units.jl") + @safetestset "LabelledArrays Test" include("labelledarrays.jl") + @safetestset "Mass Matrix Test" include("mass_matrix.jl") + @safetestset "SteadyStateSystem Test" include("steadystatesystems.jl") + @safetestset "SDESystem Test" include("sdesystem.jl") + @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl") + @safetestset "InitializationSystem Test" include("initializationsystem.jl") + @safetestset "PDE Construction Test" include("pde.jl") + @safetestset "JumpSystem Test" include("jumpsystem.jl") + @safetestset "Constraints Test" include("constraints.jl") + @safetestset "Reduction Test" include("reduction.jl") + @safetestset "Split Parameters Test" include("split_parameters.jl") + @safetestset "StaticArrays Test" include("static_arrays.jl") + @safetestset "Components Test" include("components.jl") + @safetestset "Model Parsing Test" include("model_parsing.jl") + @safetestset "print_tree" include("print_tree.jl") + @safetestset "Error Handling" include("error_handling.jl") + @safetestset "StructuralTransformations" include("structural_transformation/runtests.jl") + @safetestset "State Selection Test" include("state_selection.jl") + @safetestset "Symbolic Event Test" include("symbolic_events.jl") + @safetestset "Stream Connect Test" include("stream_connectors.jl") + @safetestset "Domain Connect Test" include("domain_connectors.jl") + @safetestset "Lowering Integration Test" include("lowering_solving.jl") + @safetestset "Test Big System Usage" include("bigsystem.jl") + @safetestset "Dependency Graph Test" include("dep_graphs.jl") + @safetestset "Function Registration Test" include("function_registration.jl") + @safetestset "Precompiled Modules Test" include("precompile_test.jl") + @safetestset "Variable Utils Test" include("variable_utils.jl") + @safetestset "Variable Metadata Test" include("test_variable_metadata.jl") + @safetestset "DAE Jacobians Test" include("dae_jacobian.jl") + @safetestset "Jacobian Sparsity" include("jacobiansparsity.jl") + @safetestset "Modelingtoolkitize Test" include("modelingtoolkitize.jl") + @safetestset "OptimizationSystem Test" include("optimizationsystem.jl") + @safetestset "FuncAffect Test" include("funcaffect.jl") + @safetestset "Constants Test" include("constants.jl") + @safetestset "Parameter Dependency Test" include("parameter_dependencies.jl") + @safetestset "Generate Custom Function Test" include("generate_custom_function.jl") + end end if GROUP == "All" || GROUP == "InterfaceII" From 3ea493aa9504aa46d6cca2f34a4dc634c3350b57 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 14:03:45 -0600 Subject: [PATCH 35/47] handle dds and t --- .../StructuralTransformations.jl | 3 +- .../symbolics_tearing.jl | 6 +++ src/systems/abstractsystem.jl | 1 + src/systems/diffeqs/abstractodesystem.jl | 38 +++++++++++++------ src/systems/diffeqs/odesystem.jl | 20 +++++++--- src/systems/nonlinear/initializesystem.jl | 23 +++++++---- src/systems/systems.jl | 8 +++- 7 files changed, 72 insertions(+), 27 deletions(-) diff --git a/src/structural_transformation/StructuralTransformations.jl b/src/structural_transformation/StructuralTransformations.jl index 5766c92e04..aa2f96da15 100644 --- a/src/structural_transformation/StructuralTransformations.jl +++ b/src/structural_transformation/StructuralTransformations.jl @@ -23,7 +23,8 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Di IncrementalCycleTracker, add_edge_checked!, topological_sort, invalidate_cache!, Substitutions, get_or_construct_tearing_state, filter_kwargs, lower_varname, setio, SparseMatrixCLIL, - fast_substitute, get_fullvars, has_equations, observed + fast_substitute, get_fullvars, has_equations, observed, + Schedule using ModelingToolkit.BipartiteGraphs import .BipartiteGraphs: invview, complete diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index 30471c7760..eb9d3ddfd6 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -555,6 +555,12 @@ function tearing_reassemble(state::TearingState, var_eq_matching; # TODO: compute the dependency correctly so that we don't have to do this obs = [fast_substitute(observed(sys), obs_sub); subeqs] @set! sys.observed = obs + + # Only makes sense for time-dependent + # TODO: generalize to SDE + if sys isa ODESystem + @set! sys.schedule = Schedule(var_eq_matching, dummy_sub) + end @set! state.sys = sys @set! sys.tearing_state = state return invalidate_cache!(sys) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 3ff6527755..a62aa4bce8 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -576,6 +576,7 @@ for prop in [:eqs :preface :torn_matching :initializesystem + :schedule :tearing_state :substitutions :metadata diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 430e9e5a85..4e999de2de 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1,3 +1,8 @@ +struct Schedule + var_eq_matching + dummy_sub +end + function filter_kwargs(kwargs) kwargs = Dict(kwargs) for key in keys(kwargs) @@ -326,7 +331,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, expression_module = eval_module, checkbounds = checkbounds, kwargs...) f_oop, f_iip = eval_expression ? - (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) : + ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) : f_gen f(u, p, t) = f_oop(u, p, t) f(du, u, p, t) = f_iip(du, u, p, t) @@ -351,7 +356,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, expression_module = eval_module, checkbounds = checkbounds, kwargs...) tgrad_oop, tgrad_iip = eval_expression ? - (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in tgrad_gen) : + ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in tgrad_gen) : tgrad_gen if p isa Tuple __tgrad(u, p, t) = tgrad_oop(u, p..., t) @@ -373,7 +378,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, expression_module = eval_module, checkbounds = checkbounds, kwargs...) jac_oop, jac_iip = eval_expression ? - (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in jac_gen) : + ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in jac_gen) : jac_gen _jac(u, p, t) = jac_oop(u, p, t) _jac(J, u, p, t) = jac_iip(J, u, p, t) @@ -541,7 +546,7 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) expression_module = eval_module, checkbounds = checkbounds, kwargs...) f_oop, f_iip = eval_expression ? - (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) : + ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) : f_gen f(du, u, p, t) = f_oop(du, u, p, t) f(du, u, p::MTKParameters, t) = f_oop(du, u, p..., t) @@ -555,7 +560,7 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) expression_module = eval_module, checkbounds = checkbounds, kwargs...) jac_oop, jac_iip = eval_expression ? - (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in jac_gen) : + ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in jac_gen) : jac_gen _jac(du, u, p, ˍ₋gamma, t) = jac_oop(du, u, p, ˍ₋gamma, t) _jac(du, u, p::MTKParameters, ˍ₋gamma, t) = jac_oop(du, u, p..., ˍ₋gamma, t) @@ -624,7 +629,7 @@ function DiffEqBase.DDEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) expression = Val{true}, expression_module = eval_module, checkbounds = checkbounds, kwargs...) - f_oop, f_iip = (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) + f_oop, f_iip = ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) f(u, h, p, t) = f_oop(u, h, p, t) f(u, h, p::MTKParameters, t) = f_oop(u, h, p..., t) f(du, u, h, p, t) = f_iip(du, u, h, p, t) @@ -649,10 +654,10 @@ function DiffEqBase.SDDEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys expression = Val{true}, expression_module = eval_module, checkbounds = checkbounds, kwargs...) - f_oop, f_iip = (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) + f_oop, f_iip = ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) g_gen = generate_diffusion_function(sys, dvs, ps; expression = Val{true}, isdde = true, kwargs...) - g_oop, g_iip = (drop_expr(@RuntimeGeneratedFunction(ex)) for ex in g_gen) + g_oop, g_iip = ((@RuntimeGeneratedFunction(ex)) for ex in g_gen) f(u, h, p, t) = f_oop(u, h, p, t) f(u, h, p::MTKParameters, t) = f_oop(u, h, p..., t) f(du, u, h, p, t) = f_iip(du, u, h, p, t) @@ -849,6 +854,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; symbolic_u0 = false, u0_constructor = identity, guesses = Dict(), + t = nothing, warn_initialize_determined = true, kwargs...) eqs = equations(sys) @@ -870,7 +876,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; u0map = unknowns(sys) .=> u0map end initializeprob = ModelingToolkit.InitializationProblem( - sys, u0map, parammap; guesses, warn_initialize_determined) + sys, t, u0map, parammap; guesses, warn_initialize_determined) initializeprobmap = getu(initializeprob, unknowns(sys)) zerovars = setdiff(unknowns(sys), keys(defaults(sys))) .=> 0.0 @@ -1101,6 +1107,7 @@ function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan end f, du0, u0, p = process_DEProblem(DAEFunction{iip}, sys, u0map, parammap; implicit_dae = true, du0map = du0map, check_length, + t = tspan !== nothing ? tspan[1] : tspan, warn_initialize_determined, kwargs...) diffvars = collect_differential_variables(sys) sts = unknowns(sys) @@ -1277,6 +1284,7 @@ function ODEProblemExpr{iip}(sys::AbstractODESystem, u0map, tspan, error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `ODEProblemExpr`") end f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap; check_length, + t = tspan !== nothing ? tspan[1] : tspan, kwargs...) linenumbers = get(kwargs, :linenumbers, true) kwargs = filter_kwargs(kwargs) @@ -1322,6 +1330,7 @@ function DAEProblemExpr{iip}(sys::AbstractODESystem, du0map, u0map, tspan, error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEProblemExpr`") end f, du0, u0, p = process_DEProblem(DAEFunctionExpr{iip}, sys, u0map, parammap; + t = tspan !== nothing ? tspan[1] : tspan, implicit_dae = true, du0map = du0map, check_length, kwargs...) linenumbers = get(kwargs, :linenumbers, true) @@ -1505,11 +1514,11 @@ function InitializationProblem(sys::AbstractODESystem, args...; kwargs...) InitializationProblem{true}(sys, args...; kwargs...) end -function InitializationProblem(sys::AbstractODESystem, +function InitializationProblem(sys::AbstractODESystem, t, u0map::StaticArray, args...; kwargs...) - InitializationProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) + InitializationProblem{false, SciMLBase.FullSpecialize}(sys, t, u0map, args...; kwargs...) end function InitializationProblem{true}(sys::AbstractODESystem, args...; kwargs...) @@ -1520,7 +1529,8 @@ function InitializationProblem{false}(sys::AbstractODESystem, args...; kwargs... InitializationProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) end -function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], +function InitializationProblem{iip, specialize}(sys::AbstractODESystem, + t, u0map = [], parammap = DiffEqBase.NullParameters(); guesses = [], check_length = true, @@ -1530,6 +1540,8 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0map = error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`") end + @show u0map + if isempty(u0map) && get_initializesystem(sys) !== nothing isys = get_initializesystem(sys) elseif isempty(u0map) && get_initializesystem(sys) === nothing @@ -1548,6 +1560,8 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0map = if warn_initialize_determined && neqs < nunknown @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false." end + + parammap isa DiffEqBase.NullParameters ? [independent_variable(sys) => t] : merge(todict(parammap), Dict(independent_variable(sys) => t)) if neqs == nunknown NonlinearProblem(isys, guesses, parammap) diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 7ee5f705f4..2076c64191 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -101,6 +101,10 @@ struct ODESystem <: AbstractODESystem """ initializesystem::Union{Nothing, NonlinearSystem} """ + The schedule for the code generation process. + """ + schedule::Any + """ Type of the system. """ connector_type::Any @@ -165,11 +169,13 @@ struct ODESystem <: AbstractODESystem """ parent::Any + function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, - torn_matching, initializesystem, connector_type, preface, cevents, - devents, parameter_dependencies, metadata = nothing, gui_metadata = nothing, - tearing_state = nothing, + torn_matching, initializesystem, schedule, connector_type, preface, cevents, + devents, parameter_dependencies, + metadata = nothing, gui_metadata = nothing, + tearing_state = nothing, substitutions = nothing, complete = false, index_cache = nothing, discrete_subsystems = nothing, solved_unknowns = nothing, split_idxs = nothing, parent = nothing; checks::Union{Bool, Int} = true) @@ -185,7 +191,8 @@ struct ODESystem <: AbstractODESystem end new(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, torn_matching, - initializesystem, connector_type, preface, cevents, devents, parameter_dependencies, metadata, + initializesystem, schedule, connector_type, preface, cevents, devents, parameter_dependencies, + metadata, gui_metadata, tearing_state, substitutions, complete, index_cache, discrete_subsystems, solved_unknowns, split_idxs, parent) end @@ -202,12 +209,13 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; defaults = _merge(Dict(default_u0), Dict(default_p)), guesses = Dict(), initializesystem = nothing, + schedule = nothing, connector_type = nothing, preface = nothing, continuous_events = nothing, discrete_events = nothing, parameter_dependencies = nothing, - checks = true, + checks = true, metadata = nothing, gui_metadata = nothing) name === nothing && @@ -253,7 +261,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; ODESystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)), deqs, iv′, dvs′, ps′, tspan, var_to_name, ctrl′, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, nothing, initializesystem, - connector_type, preface, cont_callbacks, disc_callbacks, parameter_dependencies, + schedule, connector_type, preface, cont_callbacks, disc_callbacks, parameter_dependencies, metadata, gui_metadata, checks = checks) end diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 51e9c49480..4f37f375fb 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -6,7 +6,9 @@ Generate `NonlinearSystem` which initializes an ODE problem from specified initi function generate_initializesystem(sys::ODESystem; u0map = Dict(), name = nameof(sys), - guesses = Dict(), check_defguess = false, kwargs...) + guesses = Dict(), check_defguess = false, + default_dd_value = 0.0, + kwargs...) sts, eqs = unknowns(sys), equations(sys) idxs_diff = isdiffeq.(eqs) idxs_alge = .!idxs_diff @@ -18,11 +20,18 @@ function generate_initializesystem(sys::ODESystem; defs = merge(defaults(sys), todict(u0map)) full_states = [sts; getfield.((observed(sys)), :lhs)] + set_full_states = Set(full_states) + guesses = todict(guesses) + schedule = getfield(sys, :schedule) - # Refactor to ODESystem construction - # should be ModelingToolkit.guesses(sys) + dd_guess = if schedule !== nothing + guessmap = [x[2]=>get(guesses, x[1], default_dd_value) for x in schedule.dummy_sub] + Dict(filter(x->!isnothing(x[1]) && x[1]∈set_full_states,guessmap)) + else + Dict() + end - guesses = merge(get_guesses(sys), todict(guesses)) + guesses = merge(get_guesses(sys), todict(guesses), dd_guess) for st in full_states if st ∈ keys(defs) @@ -44,13 +53,13 @@ function generate_initializesystem(sys::ODESystem; end end - pars = parameters(sys) + pars = [parameters(sys); independent_variable(sys)] nleqs = [eqs_ics; observed(sys)] - + sys_nl = NonlinearSystem(nleqs, full_states, pars; - defaults = merge(ModelingToolkit.defaults(sys), todict(u0)), + defaults = merge(ModelingToolkit.defaults(sys), todict(u0), dd_guess), name, kwargs...) diff --git a/src/systems/systems.jl b/src/systems/systems.jl index d473edc92b..9a08a89663 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -26,7 +26,13 @@ function structural_simplify( else newsys = newsys′ end - @set! newsys.parent = complete(sys; split) + if newsys isa ODESystem + schedule = newsys.schedule + @set! newsys.parent = complete(sys; split) + @set! newsys.schedule = schedule + else + @set! newsys.parent = complete(sys; split) + end newsys = complete(newsys; split) if newsys′ isa Tuple idxs = [parameter_index(newsys, i) for i in io[1]] From 734e19575680588dc8300cc37c00b635ad46f4ae Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 14:04:52 -0600 Subject: [PATCH 36/47] format --- src/structural_transformation/symbolics_tearing.jl | 2 +- src/systems/diffeqs/abstractodesystem.jl | 14 ++++++++------ src/systems/diffeqs/odesystem.jl | 7 +++---- src/systems/nonlinear/initializesystem.jl | 9 +++++---- src/systems/systems.jl | 2 +- 5 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index eb9d3ddfd6..23a5258104 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -560,7 +560,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching; # TODO: generalize to SDE if sys isa ODESystem @set! sys.schedule = Schedule(var_eq_matching, dummy_sub) - end + end @set! state.sys = sys @set! sys.tearing_state = state return invalidate_cache!(sys) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 4e999de2de..d327a96632 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1,6 +1,6 @@ struct Schedule - var_eq_matching - dummy_sub + var_eq_matching::Any + dummy_sub::Any end function filter_kwargs(kwargs) @@ -1518,7 +1518,8 @@ function InitializationProblem(sys::AbstractODESystem, t, u0map::StaticArray, args...; kwargs...) - InitializationProblem{false, SciMLBase.FullSpecialize}(sys, t, u0map, args...; kwargs...) + InitializationProblem{false, SciMLBase.FullSpecialize}( + sys, t, u0map, args...; kwargs...) end function InitializationProblem{true}(sys::AbstractODESystem, args...; kwargs...) @@ -1529,7 +1530,7 @@ function InitializationProblem{false}(sys::AbstractODESystem, args...; kwargs... InitializationProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) end -function InitializationProblem{iip, specialize}(sys::AbstractODESystem, +function InitializationProblem{iip, specialize}(sys::AbstractODESystem, t, u0map = [], parammap = DiffEqBase.NullParameters(); guesses = [], @@ -1560,8 +1561,9 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, if warn_initialize_determined && neqs < nunknown @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false." end - - parammap isa DiffEqBase.NullParameters ? [independent_variable(sys) => t] : merge(todict(parammap), Dict(independent_variable(sys) => t)) + + parammap isa DiffEqBase.NullParameters ? [independent_variable(sys) => t] : + merge(todict(parammap), Dict(independent_variable(sys) => t)) if neqs == nunknown NonlinearProblem(isys, guesses, parammap) diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 2076c64191..da973a6f46 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -169,13 +169,12 @@ struct ODESystem <: AbstractODESystem """ parent::Any - function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, torn_matching, initializesystem, schedule, connector_type, preface, cevents, devents, parameter_dependencies, metadata = nothing, gui_metadata = nothing, - tearing_state = nothing, + tearing_state = nothing, substitutions = nothing, complete = false, index_cache = nothing, discrete_subsystems = nothing, solved_unknowns = nothing, split_idxs = nothing, parent = nothing; checks::Union{Bool, Int} = true) @@ -191,7 +190,7 @@ struct ODESystem <: AbstractODESystem end new(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, torn_matching, - initializesystem, schedule, connector_type, preface, cevents, devents, parameter_dependencies, + initializesystem, schedule, connector_type, preface, cevents, devents, parameter_dependencies, metadata, gui_metadata, tearing_state, substitutions, complete, index_cache, discrete_subsystems, solved_unknowns, split_idxs, parent) @@ -215,7 +214,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; continuous_events = nothing, discrete_events = nothing, parameter_dependencies = nothing, - checks = true, + checks = true, metadata = nothing, gui_metadata = nothing) name === nothing && diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 4f37f375fb..97e5d29467 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -6,7 +6,7 @@ Generate `NonlinearSystem` which initializes an ODE problem from specified initi function generate_initializesystem(sys::ODESystem; u0map = Dict(), name = nameof(sys), - guesses = Dict(), check_defguess = false, + guesses = Dict(), check_defguess = false, default_dd_value = 0.0, kwargs...) sts, eqs = unknowns(sys), equations(sys) @@ -25,8 +25,9 @@ function generate_initializesystem(sys::ODESystem; schedule = getfield(sys, :schedule) dd_guess = if schedule !== nothing - guessmap = [x[2]=>get(guesses, x[1], default_dd_value) for x in schedule.dummy_sub] - Dict(filter(x->!isnothing(x[1]) && x[1]∈set_full_states,guessmap)) + guessmap = [x[2] => get(guesses, x[1], default_dd_value) + for x in schedule.dummy_sub] + Dict(filter(x -> !isnothing(x[1]) && x[1] ∈ set_full_states, guessmap)) else Dict() end @@ -55,7 +56,7 @@ function generate_initializesystem(sys::ODESystem; pars = [parameters(sys); independent_variable(sys)] nleqs = [eqs_ics; observed(sys)] - + sys_nl = NonlinearSystem(nleqs, full_states, pars; diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 9a08a89663..2993fe0089 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -26,7 +26,7 @@ function structural_simplify( else newsys = newsys′ end - if newsys isa ODESystem + if newsys isa ODESystem schedule = newsys.schedule @set! newsys.parent = complete(sys; split) @set! newsys.schedule = schedule From 90ceeda694924193d3d79da57bb924ca96babd2a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 15:17:55 -0600 Subject: [PATCH 37/47] update initialization problem usage for t0 --- src/systems/diffeqs/abstractodesystem.jl | 2 +- test/initializationsystem.jl | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index d327a96632..5be3079ac1 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1531,7 +1531,7 @@ function InitializationProblem{false}(sys::AbstractODESystem, args...; kwargs... end function InitializationProblem{iip, specialize}(sys::AbstractODESystem, - t, u0map = [], + t::Number, u0map = [], parammap = DiffEqBase.NullParameters(); guesses = [], check_length = true, diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index 2e6fc8c95c..157788437b 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -17,7 +17,7 @@ sol = solve(initprob) @test SciMLBase.successful_retcode(sol) @test maximum(abs.(sol[conditions])) < 1e-14 -initprob = ModelingToolkit.InitializationProblem(pend, [x => 1, y => 0], [g => 1]; +initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [x => 1, y => 0], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) @test initprob isa NonlinearProblem sol = solve(initprob) @@ -26,7 +26,7 @@ sol = solve(initprob) @test maximum(abs.(sol[conditions])) < 1e-14 initprob = ModelingToolkit.InitializationProblem( - pend, [], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) + pend, 0.0, [], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) @test initprob isa NonlinearLeastSquaresProblem sol = solve(initprob) @test !SciMLBase.successful_retcode(sol) @@ -216,7 +216,7 @@ end end @mtkbuild sys = System() -initprob = ModelingToolkit.InitializationProblem(sys) +initprob = ModelingToolkit.InitializationProblem(sys, 0.0) conditions = getfield.(equations(initprob.f.sys), :rhs) @test initprob isa NonlinearLeastSquaresProblem @@ -296,7 +296,7 @@ end end @mtkbuild sys = MassDamperSystem() -initprob = ModelingToolkit.InitializationProblem(sys) +initprob = ModelingToolkit.InitializationProblem(sys, 0.0) @test initprob isa NonlinearProblem initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12) @test SciMLBase.successful_retcode(initsol) From 920fcc2151ead7a190392a497050002e3bb440af Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 15:27:02 -0600 Subject: [PATCH 38/47] Remove early caching initialization system --- src/systems/diffeqs/abstractodesystem.jl | 2 -- src/systems/systemstructure.jl | 18 ------------------ 2 files changed, 20 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 5be3079ac1..77f627641f 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1541,8 +1541,6 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`") end - @show u0map - if isempty(u0map) && get_initializesystem(sys) !== nothing isys = get_initializesystem(sys) elseif isempty(u0map) && get_initializesystem(sys) === nothing diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 3f022640fe..6c94751b5d 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -636,23 +636,5 @@ function _structural_simplify!(state::TearingState, io; simplify = false, fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)] @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns) - ci = infer_clocks!(ClockInference(state)) - # TODO: make it work with clocks - if sys isa ODESystem && all(isequal(Continuous()), ci.var_domain) && - (!has_io || - !all(all(x -> !(typeof(x) <: Union{Sample, Hold, ShiftIndex}), io))) - isys = ModelingToolkit.generate_initializesystem(sys) - !isempty(equations(isys)) && - (isys = structural_simplify(isys; fully_determined = false)) - @set! sys.initializesystem = isys - neqs = length(equations(isys)) - nunknown = length(unknowns(isys)) - if warn_initialize_determined && neqs > nunknown - @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares if $(nunknown - neqs) defaults are not supplied at construction time." - end - if warn_initialize_determined && neqs < nunknown - @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares" - end - end ModelingToolkit.invalidate_cache!(sys), input_idxs end From 0b6edee752a51346133b7feb10eeae088a0cf4ec Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 15:51:34 -0600 Subject: [PATCH 39/47] guesses --- src/systems/nonlinear/initializesystem.jl | 2 +- test/state_selection.jl | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 97e5d29467..d486c4922c 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -27,7 +27,7 @@ function generate_initializesystem(sys::ODESystem; dd_guess = if schedule !== nothing guessmap = [x[2] => get(guesses, x[1], default_dd_value) for x in schedule.dummy_sub] - Dict(filter(x -> !isnothing(x[1]) && x[1] ∈ set_full_states, guessmap)) + Dict(filter(x -> !isnothing(x[1]), guessmap)) else Dict() end diff --git a/test/state_selection.jl b/test/state_selection.jl index c1e4328715..e706c48d2d 100644 --- a/test/state_selection.jl +++ b/test/state_selection.jl @@ -185,9 +185,10 @@ let mo_1 => 0 mo_2 => 1 mo_3 => 2 + Ek_2 => 2 Ek_3 => 3] - prob1 = ODEProblem(sys, u0, (0.0, 0.1)) - prob2 = ODEProblem(sys, u0, (0.0, 0.1)) + prob1 = ODEProblem(sys, [], (0.0, 0.1), guesses = u0) + prob2 = ODEProblem(sys, [], (0.0, 0.1), guesses = u0) @test solve(prob1, FBDF()).retcode == ReturnCode.Success @test solve(prob2, FBDF()).retcode == ReturnCode.Success end From d78c68011f16e5da9a40821b1b895e0134c82133 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 15:58:03 -0600 Subject: [PATCH 40/47] handle empty parammap --- src/systems/diffeqs/abstractodesystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 77f627641f..a07deb76b1 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1560,7 +1560,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false." end - parammap isa DiffEqBase.NullParameters ? [independent_variable(sys) => t] : + parammap isa DiffEqBase.NullParameters || isempty(parammap) ? [independent_variable(sys) => t] : merge(todict(parammap), Dict(independent_variable(sys) => t)) if neqs == nunknown From 79eca95d3428a9ad07f03d3949dae9fc0246704b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 15:59:41 -0600 Subject: [PATCH 41/47] format --- src/systems/diffeqs/abstractodesystem.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index a07deb76b1..e403f47ceb 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1560,7 +1560,8 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false." end - parammap isa DiffEqBase.NullParameters || isempty(parammap) ? [independent_variable(sys) => t] : + parammap isa DiffEqBase.NullParameters || isempty(parammap) ? + [independent_variable(sys) => t] : merge(todict(parammap), Dict(independent_variable(sys) => t)) if neqs == nunknown From ab935d73aaa6887018966324110822ef12b858e7 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 16:11:10 -0600 Subject: [PATCH 42/47] Fix schedule setting --- src/systems/systems.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 2993fe0089..3987ae89ed 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -27,9 +27,7 @@ function structural_simplify( newsys = newsys′ end if newsys isa ODESystem - schedule = newsys.schedule @set! newsys.parent = complete(sys; split) - @set! newsys.schedule = schedule else @set! newsys.parent = complete(sys; split) end From 2ccf8fa9bd94422981fcb9427e387b1100fb0424 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 27 Feb 2024 21:44:54 -0600 Subject: [PATCH 43/47] a few fixes --- src/systems/nonlinear/initializesystem.jl | 2 +- test/initializationsystem.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index d486c4922c..827d5ca2e1 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -54,7 +54,7 @@ function generate_initializesystem(sys::ODESystem; end end - pars = [parameters(sys); independent_variable(sys)] + pars = [parameters(sys); get_iv(sys)] nleqs = [eqs_ics; observed(sys)] sys_nl = NonlinearSystem(nleqs, diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index 157788437b..69032eeb3b 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -8,7 +8,7 @@ eqs = [D(D(x)) ~ λ * x x^2 + y^2 ~ 1] @mtkbuild pend = ODESystem(eqs, t) -initprob = ModelingToolkit.InitializationProblem(pend, [], [g => 1]; +initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [], [g => 1]; guesses = [ModelingToolkit.missing_variable_defaults(pend); x => 1; y => 0.2]) conditions = getfield.(equations(initprob.f.sys), :rhs) From 9cf1768d53f7b8c6fd8f2aab9c64361c86be917b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 28 Feb 2024 07:18:17 -0600 Subject: [PATCH 44/47] Fix test cases --- test/components.jl | 6 ++--- test/state_selection.jl | 24 +++++++++---------- .../index_reduction.jl | 3 +-- 3 files changed, 16 insertions(+), 17 deletions(-) diff --git a/test/components.jl b/test/components.jl index 233d1c78c5..a7d69ab24b 100644 --- a/test/components.jl +++ b/test/components.jl @@ -98,9 +98,9 @@ let @named rc_model2 = compose(_rc_model2, [resistor, resistor2, capacitor, source, ground]) sys2 = structural_simplify(rc_model2) - prob2 = ODEProblem(sys2, u0, (0, 10.0)) + prob2 = ODEProblem(sys2, [], (0, 10.0), guesses = u0) sol2 = solve(prob2, Rosenbrock23()) - @test sol2[source.p.i] ≈ sol2[rc_model2.source.p.i] ≈ -sol2[capacitor.i] + @test sol2[source.p.i] ≈ sol2[rc_model2.source.p.i] ≈ sol2[capacitor.i] end # Outer/inner connections @@ -155,7 +155,7 @@ include("../examples/serial_inductor.jl") sys = structural_simplify(ll_model) @test length(equations(sys)) == 2 u0 = unknowns(sys) .=> 0 -@test_nowarn ODEProblem(sys, u0, (0, 10.0)) +@test_nowarn ODEProblem(sys, [], (0, 10.0), guesses = u0, warn_initialize_determined = false) prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0, u0, (0, 0.5)) sol = solve(prob, DFBDF()) @test sol.retcode == SciMLBase.ReturnCode.Success diff --git a/test/state_selection.jl b/test/state_selection.jl index e706c48d2d..8a198fc2ef 100644 --- a/test/state_selection.jl +++ b/test/state_selection.jl @@ -45,16 +45,16 @@ end # 1516 let @connector function Fluid_port(; name, p = 101325.0, m = 0.0, T = 293.15) - sts = @variables p(t)=p m(t)=m [connect = Flow] T(t)=T [connect = Stream] + sts = @variables p(t) [guess=p] m(t) [guess = m, connect = Flow] T(t) [guess=T, connect = Stream] ODESystem(Equation[], t, sts, []; name = name) end - + #this one is for latter @connector function Heat_port(; name, Q = 0.0, T = 293.15) - sts = @variables T(t)=T Q(t)=Q [connect = Flow] + sts = @variables T(t) [guess=T] Q(t) [guess = Q, connect = Flow] ODESystem(Equation[], t, sts, []; name = name) end - + # like ground but for fluid systems (fluid_port.m is expected to be zero in closed loop) function Compensator(; name, p = 101325.0, T_back = 273.15) @named fluid_port = Fluid_port() @@ -63,7 +63,7 @@ let fluid_port.T ~ T_back] compose(ODESystem(eqs, t, [], ps; name = name), fluid_port) end - + function Source(; name, delta_p = 100, T_feed = 293.15) @named supply_port = Fluid_port() # expected to feed connected pipe -> m<0 @named return_port = Fluid_port() # expected to receive from connected pipe -> m>0 @@ -74,7 +74,7 @@ let return_port.T ~ T_feed] compose(ODESystem(eqs, t, [], ps; name = name), [supply_port, return_port]) end - + function Substation(; name, T_return = 343.15) @named supply_port = Fluid_port() # expected to receive from connected pipe -> m>0 @named return_port = Fluid_port() # expected to feed connected pipe -> m<0 @@ -85,12 +85,12 @@ let return_port.T ~ T_return] compose(ODESystem(eqs, t, [], ps; name = name), [supply_port, return_port]) end - + function Pipe(; name, L = 1000, d = 0.1, N = 100, rho = 1000, f = 1) @named fluid_port_a = Fluid_port() @named fluid_port_b = Fluid_port() ps = @parameters L=L d=d rho=rho f=f N=N - sts = @variables v(t)=0.0 dp_z(t)=0.0 + sts = @variables v(t) [guess=0.0] dp_z(t) [guess=0.0] eqs = [fluid_port_a.m ~ -fluid_port_b.m fluid_port_a.T ~ instream(fluid_port_a.T) fluid_port_b.T ~ fluid_port_a.T @@ -114,7 +114,7 @@ let connect(return_pipe.fluid_port_a, source.return_port)] compose(ODESystem(eqs, t, [], ps; name = name), subs) end - + @named system = System(L = 10) @unpack supply_pipe, return_pipe = system sys = structural_simplify(system) @@ -122,10 +122,10 @@ let system.supply_pipe.v => 0.1, system.return_pipe.v => 0.1, D(supply_pipe.v) => 0.0, 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 = ODEProblem(sys, u0, (0.0, 10.0), []) + prob1 = ODEProblem(sys, [], (0.0, 10.0), [], guesses = u0) + prob2 = ODEProblem(sys, [], (0.0, 10.0), [], guesses = u0) prob3 = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, u0, (0.0, 10.0), []) - @test solve(prob1, FBDF()).retcode == ReturnCode.Success + @test solve(prob1, FBDF()).retcode == ReturnCode.Success #@test solve(prob2, FBDF()).retcode == ReturnCode.Success @test solve(prob3, DFBDF()).retcode == ReturnCode.Success end diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl index 8dc06d680c..7b00c1c317 100644 --- a/test/structural_transformation/index_reduction.jl +++ b/test/structural_transformation/index_reduction.jl @@ -144,14 +144,13 @@ let sys = structural_simplify(pendulum2) D(D(y)) => 0.0, x => sqrt(2) / 2, y => sqrt(2) / 2, - T => 0.0 ] p = [ L => 1.0, g => 9.8 ] - prob_auto = ODEProblem(sys, u0, (0.0, 0.5), p) + prob_auto = ODEProblem(sys, u0, (0.0, 0.5), p, guesses = [T => 0.0]) sol = solve(prob_auto, FBDF()) @test sol.retcode == ReturnCode.Success @test norm(sol[x] .^ 2 + sol[y] .^ 2 .- 1) < 1e-2 From 537fb7c6ae537c54e2fe4a0b971984b00432b115 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 28 Feb 2024 07:21:20 -0600 Subject: [PATCH 45/47] format --- test/components.jl | 3 ++- test/state_selection.jl | 21 ++++++++++--------- .../index_reduction.jl | 2 +- 3 files changed, 14 insertions(+), 12 deletions(-) diff --git a/test/components.jl b/test/components.jl index a7d69ab24b..76f456887e 100644 --- a/test/components.jl +++ b/test/components.jl @@ -155,7 +155,8 @@ include("../examples/serial_inductor.jl") sys = structural_simplify(ll_model) @test length(equations(sys)) == 2 u0 = unknowns(sys) .=> 0 -@test_nowarn ODEProblem(sys, [], (0, 10.0), guesses = u0, warn_initialize_determined = false) +@test_nowarn ODEProblem( + sys, [], (0, 10.0), guesses = u0, warn_initialize_determined = false) prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0, u0, (0, 0.5)) sol = solve(prob, DFBDF()) @test sol.retcode == SciMLBase.ReturnCode.Success diff --git a/test/state_selection.jl b/test/state_selection.jl index 8a198fc2ef..c847955a85 100644 --- a/test/state_selection.jl +++ b/test/state_selection.jl @@ -45,16 +45,17 @@ end # 1516 let @connector function Fluid_port(; name, p = 101325.0, m = 0.0, T = 293.15) - sts = @variables p(t) [guess=p] m(t) [guess = m, connect = Flow] T(t) [guess=T, connect = Stream] + sts = @variables p(t) [guess = p] m(t) [guess = m, connect = Flow] T(t) [ + guess = T, connect = Stream] ODESystem(Equation[], t, sts, []; name = name) end - + #this one is for latter @connector function Heat_port(; name, Q = 0.0, T = 293.15) - sts = @variables T(t) [guess=T] Q(t) [guess = Q, connect = Flow] + sts = @variables T(t) [guess = T] Q(t) [guess = Q, connect = Flow] ODESystem(Equation[], t, sts, []; name = name) end - + # like ground but for fluid systems (fluid_port.m is expected to be zero in closed loop) function Compensator(; name, p = 101325.0, T_back = 273.15) @named fluid_port = Fluid_port() @@ -63,7 +64,7 @@ let fluid_port.T ~ T_back] compose(ODESystem(eqs, t, [], ps; name = name), fluid_port) end - + function Source(; name, delta_p = 100, T_feed = 293.15) @named supply_port = Fluid_port() # expected to feed connected pipe -> m<0 @named return_port = Fluid_port() # expected to receive from connected pipe -> m>0 @@ -74,7 +75,7 @@ let return_port.T ~ T_feed] compose(ODESystem(eqs, t, [], ps; name = name), [supply_port, return_port]) end - + function Substation(; name, T_return = 343.15) @named supply_port = Fluid_port() # expected to receive from connected pipe -> m>0 @named return_port = Fluid_port() # expected to feed connected pipe -> m<0 @@ -85,12 +86,12 @@ let return_port.T ~ T_return] compose(ODESystem(eqs, t, [], ps; name = name), [supply_port, return_port]) end - + function Pipe(; name, L = 1000, d = 0.1, N = 100, rho = 1000, f = 1) @named fluid_port_a = Fluid_port() @named fluid_port_b = Fluid_port() ps = @parameters L=L d=d rho=rho f=f N=N - sts = @variables v(t) [guess=0.0] dp_z(t) [guess=0.0] + sts = @variables v(t) [guess = 0.0] dp_z(t) [guess = 0.0] eqs = [fluid_port_a.m ~ -fluid_port_b.m fluid_port_a.T ~ instream(fluid_port_a.T) fluid_port_b.T ~ fluid_port_a.T @@ -114,7 +115,7 @@ let connect(return_pipe.fluid_port_a, source.return_port)] compose(ODESystem(eqs, t, [], ps; name = name), subs) end - + @named system = System(L = 10) @unpack supply_pipe, return_pipe = system sys = structural_simplify(system) @@ -125,7 +126,7 @@ let prob1 = ODEProblem(sys, [], (0.0, 10.0), [], guesses = u0) prob2 = ODEProblem(sys, [], (0.0, 10.0), [], guesses = u0) prob3 = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, u0, (0.0, 10.0), []) - @test solve(prob1, FBDF()).retcode == ReturnCode.Success + @test solve(prob1, FBDF()).retcode == ReturnCode.Success #@test solve(prob2, FBDF()).retcode == ReturnCode.Success @test solve(prob3, DFBDF()).retcode == ReturnCode.Success end diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl index 7b00c1c317..6d7c153775 100644 --- a/test/structural_transformation/index_reduction.jl +++ b/test/structural_transformation/index_reduction.jl @@ -143,7 +143,7 @@ let sys = structural_simplify(pendulum2) D(y) => 0.0, D(D(y)) => 0.0, x => sqrt(2) / 2, - y => sqrt(2) / 2, + y => sqrt(2) / 2 ] p = [ L => 1.0, From 56421fb9583804216a48f4a9bc434633eee26c10 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 28 Feb 2024 07:37:25 -0600 Subject: [PATCH 46/47] return drop_expr --- src/systems/diffeqs/abstractodesystem.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index e403f47ceb..6c64f14312 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -331,7 +331,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, expression_module = eval_module, checkbounds = checkbounds, kwargs...) f_oop, f_iip = eval_expression ? - ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) : + (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) : f_gen f(u, p, t) = f_oop(u, p, t) f(du, u, p, t) = f_iip(du, u, p, t) @@ -356,7 +356,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, expression_module = eval_module, checkbounds = checkbounds, kwargs...) tgrad_oop, tgrad_iip = eval_expression ? - ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in tgrad_gen) : + (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in tgrad_gen) : tgrad_gen if p isa Tuple __tgrad(u, p, t) = tgrad_oop(u, p..., t) @@ -378,7 +378,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, expression_module = eval_module, checkbounds = checkbounds, kwargs...) jac_oop, jac_iip = eval_expression ? - ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in jac_gen) : + (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in jac_gen) : jac_gen _jac(u, p, t) = jac_oop(u, p, t) _jac(J, u, p, t) = jac_iip(J, u, p, t) @@ -546,7 +546,7 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) expression_module = eval_module, checkbounds = checkbounds, kwargs...) f_oop, f_iip = eval_expression ? - ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) : + (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) : f_gen f(du, u, p, t) = f_oop(du, u, p, t) f(du, u, p::MTKParameters, t) = f_oop(du, u, p..., t) @@ -560,7 +560,7 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) expression_module = eval_module, checkbounds = checkbounds, kwargs...) jac_oop, jac_iip = eval_expression ? - ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in jac_gen) : + (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in jac_gen) : jac_gen _jac(du, u, p, ˍ₋gamma, t) = jac_oop(du, u, p, ˍ₋gamma, t) _jac(du, u, p::MTKParameters, ˍ₋gamma, t) = jac_oop(du, u, p..., ˍ₋gamma, t) @@ -629,7 +629,7 @@ function DiffEqBase.DDEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) expression = Val{true}, expression_module = eval_module, checkbounds = checkbounds, kwargs...) - f_oop, f_iip = ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) + f_oop, f_iip = (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) f(u, h, p, t) = f_oop(u, h, p, t) f(u, h, p::MTKParameters, t) = f_oop(u, h, p..., t) f(du, u, h, p, t) = f_iip(du, u, h, p, t) @@ -654,10 +654,10 @@ function DiffEqBase.SDDEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys expression = Val{true}, expression_module = eval_module, checkbounds = checkbounds, kwargs...) - f_oop, f_iip = ((@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) + f_oop, f_iip = (drop_expr(@RuntimeGeneratedFunction(eval_module, ex)) for ex in f_gen) g_gen = generate_diffusion_function(sys, dvs, ps; expression = Val{true}, isdde = true, kwargs...) - g_oop, g_iip = ((@RuntimeGeneratedFunction(ex)) for ex in g_gen) + g_oop, g_iip = (drop_expr(@RuntimeGeneratedFunction(ex)) for ex in g_gen) f(u, h, p, t) = f_oop(u, h, p, t) f(u, h, p::MTKParameters, t) = f_oop(u, h, p..., t) f(du, u, h, p, t) = f_iip(du, u, h, p, t) From c78321fe3159fd2d046e1a3365ceba125d2ef781 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 28 Feb 2024 08:06:31 -0600 Subject: [PATCH 47/47] Fix warning --- 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 6c64f14312..3ad6bff0ff 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1561,8 +1561,8 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, end parammap isa DiffEqBase.NullParameters || isempty(parammap) ? - [independent_variable(sys) => t] : - merge(todict(parammap), Dict(independent_variable(sys) => t)) + [get_iv(sys) => t] : + merge(todict(parammap), Dict(get_iv(sys) => t)) if neqs == nunknown NonlinearProblem(isys, guesses, parammap)