diff --git a/Project.toml b/Project.toml index 716b58eb03..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.0.1" +SciMLBase = "2.28.0" SciMLStructures = "1.0" Serialization = "1" Setfield = "0.7, 0.8, 1" 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/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/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..23a5258104 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 f95cb1fe6f..a62aa4bce8 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, :initializesystem) && get_initializesystem(sys) !== nothing + @set! sys.initializesystem = complete(get_initializesystem(sys); split) + end isdefined(sys, :complete) ? (@set! sys.complete = true) : sys end @@ -551,6 +554,7 @@ for prop in [:eqs :var_to_name :ctrls :defaults + :guesses :observed :tgrad :jac @@ -571,6 +575,8 @@ for prop in [:eqs :connections :preface :torn_matching + :initializesystem + :schedule :tearing_state :substitutions :metadata @@ -933,6 +939,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(_) = [] @@ -2259,14 +2269,15 @@ 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 b6e8c50563..3ad6bff0ff 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1,3 +1,8 @@ +struct Schedule + var_eq_matching::Any + dummy_sub::Any +end + function filter_kwargs(kwargs) kwargs = Dict(kwargs) for key in keys(kwargs) @@ -316,6 +321,8 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, sparsity = false, analytic = 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 +494,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 +503,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 +535,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 +608,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...) @@ -839,18 +853,46 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; tofloat = true, symbolic_u0 = false, u0_constructor = identity, + guesses = Dict(), + t = nothing, + warn_initialize_determined = true, kwargs...) eqs = equations(sys) dvs = unknowns(sys) ps = full_parameters(sys) iv = get_iv(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! + 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) && + ModelingToolkit.get_tearing_state(sys) !== nothing + if eltype(u0map) <: Number + u0map = unknowns(sys) .=> u0map + end + initializeprob = ModelingToolkit.InitializationProblem( + sys, t, u0map, parammap; guesses, warn_initialize_determined) + initializeprobmap = getu(initializeprob, unknowns(sys)) + + zerovars = setdiff(unknowns(sys), keys(defaults(sys))) .=> 0.0 + trueinit = identity.([zerovars; u0map]) + else + initializeprob = nothing + initializeprobmap = nothing + 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, @@ -881,6 +923,8 @@ 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, + initializeprobmap = initializeprobmap, kwargs...) implicit_dae ? (f, du0, u0, p) : (f, u0, p) end @@ -984,13 +1028,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 @@ -1055,13 +1100,15 @@ 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...) + t = tspan !== nothing ? tspan[1] : tspan, + warn_initialize_determined, kwargs...) diffvars = collect_differential_variables(sys) sts = unknowns(sys) differential_vars = map(Base.Fix2(in, diffvars), sts) @@ -1237,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) @@ -1282,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) @@ -1442,3 +1491,82 @@ 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, t, + u0map::StaticArray, + args...; + kwargs...) + InitializationProblem{false, SciMLBase.FullSpecialize}( + sys, t, 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, + t::Number, 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 + + 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 + + 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 + + parammap isa DiffEqBase.NullParameters || isempty(parammap) ? + [get_iv(sys) => t] : + merge(todict(parammap), Dict(get_iv(sys) => t)) + + if neqs == nunknown + NonlinearProblem(isys, guesses, parammap) + else + NonlinearLeastSquaresProblem(isys, guesses, parammap) + end +end diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index fe52c032e0..da973a6f46 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -88,10 +88,23 @@ 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} """ + The system for performing the initialization. + """ + initializesystem::Union{Nothing, NonlinearSystem} + """ + The schedule for the code generation process. + """ + schedule::Any + """ Type of the system. """ connector_type::Any @@ -157,9 +170,10 @@ 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, - torn_matching, connector_type, preface, cevents, - devents, parameter_dependencies, metadata = nothing, gui_metadata = nothing, + 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, substitutions = nothing, complete = false, index_cache = nothing, discrete_subsystems = nothing, solved_unknowns = nothing, @@ -175,8 +189,9 @@ 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, - connector_type, preface, cevents, devents, parameter_dependencies, metadata, + ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, torn_matching, + 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 @@ -191,6 +206,9 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), + guesses = Dict(), + initializesystem = nothing, + schedule = nothing, connector_type = nothing, preface = nothing, continuous_events = nothing, @@ -217,6 +235,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,8 +259,8 @@ 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, - connector_type, preface, cont_callbacks, disc_callbacks, parameter_dependencies, + ctrl_jac, Wfact, Wfact_t, name, systems, defaults, guesses, nothing, initializesystem, + 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 49acdf0127..827d5ca2e1 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -3,51 +3,64 @@ $(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 generate_initializesystem(sys::ODESystem; + u0map = Dict(), + name = nameof(sys), + 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 + num_alge = sum(idxs_alge) - # Algebraic equations and initial guesses are unchanged - eqs_ics = similar(eqs) - u0 = Vector{Any}(undef, length(sts)) + # Start the equations list with algebraic equations + eqs_ics = eqs[idxs_alge] + u0 = Vector{Pair}(undef, 0) + defs = merge(defaults(sys), todict(u0map)) - eqs_ics[idxs_alge] .= eqs[idxs_alge] - u0[idxs_alge] .= getmetadata.(unwrap.(sts[idxs_alge]), - Symbolics.VariableDefaultValue, - nothing) + full_states = [sts; getfield.((observed(sys)), :lhs)] + set_full_states = Set(full_states) + guesses = todict(guesses) + schedule = getfield(sys, :schedule) - for idx in findall(idxs_diff) - st = sts[idx] - if !hasdefault(st) - error("Invalid setup: unknown $(st) has no default value or equation.") - end + 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]), guessmap)) + else + Dict() + end - def = getdefault(st) - if def isa Equation - if !hasguess(st) - error("Invalid setup: unknown $(st) has an initial condition equation with no guess.") - end - guess = getguess(st) - eqs_ics[idx] = def + guesses = merge(get_guesses(sys), todict(guesses), dd_guess) - u0[idx] = guess - else - eqs_ics[idx] = st ~ def + for st in full_states + if st ∈ keys(defs) + def = defs[st] - u0[idx] = def + if def isa Equation + 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]) + else + push!(eqs_ics, st ~ def) + push!(u0, st => def) + end + elseif st ∈ keys(guesses) + push!(u0, st => guesses[st]) + elseif check_defguess + error("Invalid setup: unknown $(st) has no default value or initial guess") end end - pars = parameters(sys) - sys_nl = NonlinearSystem(eqs_ics, - sts, + pars = [parameters(sys); get_iv(sys)] + nleqs = [eqs_ics; observed(sys)] + + sys_nl = NonlinearSystem(nleqs, + full_states, pars; - defaults = Dict(sts .=> u0), + defaults = merge(ModelingToolkit.defaults(sys), todict(u0), dd_guess), name, kwargs...) diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index b00efde2ab..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) @@ -289,6 +288,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))), jac_prototype = sparse ? similar(calculate_jacobian(sys, sparse = sparse), Float64) : nothing, @@ -333,12 +334,14 @@ 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 @@ -358,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...) @@ -399,6 +404,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 +473,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) diff --git a/src/systems/systems.jl b/src/systems/systems.jl index d473edc92b..3987ae89ed 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -26,7 +26,11 @@ function structural_simplify( else newsys = newsys′ end - @set! newsys.parent = complete(sys; split) + if newsys isa ODESystem + @set! newsys.parent = complete(sys; split) + 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]] diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 064ea1314d..6c94751b5d 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 = false, kwargs...) check_consistency &= fully_determined has_io = io !== nothing @@ -635,5 +635,6 @@ 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) + ModelingToolkit.invalidate_cache!(sys), input_idxs end diff --git a/test/components.jl b/test/components.jl index 233d1c78c5..76f456887e 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,8 @@ 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/initializationsystem.jl b/test/initializationsystem.jl new file mode 100644 index 0000000000..69032eeb3b --- /dev/null +++ b/test/initializationsystem.jl @@ -0,0 +1,334 @@ +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, 0.0, [], [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, 0.0, [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, 0.0, [], [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] +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, 0.0) +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(), initializealg = BrownFullBasicInit()) +# If initialized incorrectly, then it would be InitialFailure +@test sol.retcode == SciMLBase.ReturnCode.Unstable +@test maximum(abs.(initsol[conditions][1])) < 1e-14 + +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] + 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, 0.0) +@test initprob isa NonlinearProblem +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.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) 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 diff --git a/test/reduction.jl b/test/reduction.jl index 7a8979ba23..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(eqs, xs, []) +eqs = xs ~ A * xs +@named sys′ = NonlinearSystem(eqs, [xs], []) sys = structural_simplify(sys′) # issue 958 diff --git a/test/runtests.jl b/test/runtests.jl index 96bdbbf06e..371834f351 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,54 +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 "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" diff --git a/test/state_selection.jl b/test/state_selection.jl index c1e4328715..c847955a85 100644 --- a/test/state_selection.jl +++ b/test/state_selection.jl @@ -45,13 +45,14 @@ 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 @@ -90,7 +91,7 @@ let @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 @@ -122,8 +123,8 @@ 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(prob2, FBDF()).retcode == ReturnCode.Success @@ -185,9 +186,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 diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl index 8dc06d680c..6d7c153775 100644 --- a/test/structural_transformation/index_reduction.jl +++ b/test/structural_transformation/index_reduction.jl @@ -143,15 +143,14 @@ let sys = structural_simplify(pendulum2) D(y) => 0.0, D(D(y)) => 0.0, x => sqrt(2) / 2, - y => sqrt(2) / 2, - T => 0.0 + y => sqrt(2) / 2 ] 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