Skip to content

Commit

Permalink
format
Browse files Browse the repository at this point in the history
  • Loading branch information
ChrisRackauckas committed Feb 26, 2024
1 parent 67f6f24 commit 7c0c423
Show file tree
Hide file tree
Showing 4 changed files with 33 additions and 27 deletions.
3 changes: 2 additions & 1 deletion src/systems/abstractsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
20 changes: 11 additions & 9 deletions src/systems/diffeqs/abstractodesystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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))
Expand Down
4 changes: 2 additions & 2 deletions src/systems/nonlinear/initializesystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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...)
Expand All @@ -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)]

Expand Down
33 changes: 18 additions & 15 deletions test/initializationsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
@test_throws ArgumentError prob=ODEProblem(sys, u0, tspan, p, jac = true)

0 comments on commit 7c0c423

Please sign in to comment.