ModelingToolkit.jl icon indicating copy to clipboard operation
ModelingToolkit.jl copied to clipboard

ModelingToolkit generated functions don't check input dimensions

Open mzhen77 opened this issue 2 years ago • 9 comments

Code to reproduce:

using ModelingToolkit, Plots, DifferentialEquations
using ModelingToolkitStandardLibrary.Electrical

@variables t

@named resistor = Resistor(R=1.)
@named capacitor = Capacitor(C=1.)
@named source = ConstantVoltage(V=1.)
@named ground = Ground()

rc_eqs = [
    connect(source.p, resistor.p)
    connect(resistor.n, capacitor.p)
    connect(capacitor.n, source.n)
    connect(capacitor.n, ground.g)
]

@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model,
    [resistor, capacitor, source, ground])
sys = structural_simplify(rc_model)
u0 = [
    capacitor.v => 0.0
]

params = [resistor.R => 1.0, source.V => 1.0, capacitor.C => 1.0]
prob = ODAEProblem(sys, u0, (0., 10.), params)

The next line (using only 2 parameters instead of 3) gives the wrong result without any warnings: sol = solve(prob1, Tsit5(), p=[1., 1.])

julia> sol.u
9-element Vector{Vector{Float64}}:
 [0.0]
 [5.4344286e-317]
 [5.97784145e-316]
 [6.03188388e-315]
 [6.034300412e-314]
 [6.00477522454e-313]
 [5.714879592387e-312]
 [3.644387262478e-311]
 [-7.440771759894994e-309]

The correct result is:

sol = solve(prob1, Tsit5(), p=[1., 1., 1.])
julia> sol.u
19-element Vector{Vector{Float64}}:
 [0.0]
 [9.999500016666247e-5]
 [0.001099395221772342]
 [0.011038622307372232]
 [0.0738713206981718]
 ⋮
 [0.9953772824474669]
 [0.9986536797747498]
 [0.9996930416967437]
 [0.9999471971800106]
 [0.9999499280997356]

My versions:

julia> VERSION
v"1.7.2"
pkg> st
  [0c46a032] DifferentialEquations v7.1.0
  [961ee093] ModelingToolkit v8.14.0
  [16a59e39] ModelingToolkitStandardLibrary v1.3.1
  [91a5bcdd] Plots v1.29.0

mzhen77 avatar Jun 07 '22 23:06 mzhen77

@YingboMa how is it able to run with the wrong number of params? Is that inbounds stuff?

ChrisRackauckas avatar Jun 09 '22 08:06 ChrisRackauckas

@shashi can we make build_function check the dimension of the inputs even with checkbounds=false?

YingboMa avatar Jul 01 '22 05:07 YingboMa

Well this would also not allow longer inputs that used to be allowed.

shashi avatar Jul 25 '22 19:07 shashi

Should we also check the output dimension?

shashi avatar Jul 25 '22 19:07 shashi

Where is that allowed? It's not on purpose.

ChrisRackauckas avatar Jul 25 '22 19:07 ChrisRackauckas

If this is added to the generated functions can their be a kwarg to disable it in build_function? I don’t think one would want this for a one line function that is called a ton in, for example, a jump solver.

isaacsas avatar Jul 25 '22 23:07 isaacsas

Yes, the keyword arguments are splatted onwards.

ChrisRackauckas avatar Jul 26 '22 02:07 ChrisRackauckas

The default in the jump stuff should probably disable it because yes that would have terrible performance.

ChrisRackauckas avatar Jul 26 '22 02:07 ChrisRackauckas

Alright that sounds good. So by default it will be turned on, and the downstream packages can turn it off.

shashi avatar Jul 26 '22 17:07 shashi