ModellingToolkit and arrays

Specific Domains Modelling & Simulations
, ,
1

I have the following code:

# Example one: Falling mass.
using ModelingToolkit, OrdinaryDiffEq

G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/sĀ²]

# definiting the model
@independent_variables t
@variables   pos(t)[1:3]=[0.0, 0.0,  0.0]
@variables   vel(t)[1:3]=[0.0, 0.0, 50.0] 
@variables   acc(t)[1:3]=[0.0, 0.0, -9.81] 
D = Differential(t)

pos = collect(pos)
vel = collect(vel)
acc = collect(acc)

eqs = vcat(D.(pos) .~ vel,
           D.(vel) .~ acc,
           acc    .~ G_EARTH)

@named sys = ODESystem(eqs, t)
simple_sys = structural_simplify(sys)

# running the simulation
duration = 10.0
dt = 0.02
tol = 1e-6
ts    = 0:dt:duration

prob = ODEProblem(simple_sys, nothing, (0.0, duration))
@time sol = solve(prob, Rodas5(), dt=dt, abstol=tol, reltol=tol, saveat=ts)

In older versions of ModelingToolkit the lines:

pos = collect(pos)
vel = collect(vel)
acc = collect(acc)

where not needed, but now they are. Without them I get the error message:

ERROR: LoadError: type Term has no field lhs.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] throw_no_field(::Val{:Term}, s::Symbol)
   @ Unityper ~/.julia/packages/Unityper/BiodR/src/compactify.jl:505
 [3] getproperty(x::SymbolicUtils.BasicSymbolic{Any}, s::Symbol)
   @ SymbolicUtils ~/.julia/packages/Unityper/BiodR/src/compactify.jl:310
 [4] ODESystem(eqs::Vector{Any}, iv::Num; kwargs::@Kwargs{name::Symbol})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/jEIdO/src/systems/diffeqs/odesystem.jl:299
 [5] top-level scope
   @ ~/.julia/packages/ModelingToolkit/jEIdO/src/systems/abstractsystem.jl:1810
 [6] include(fname::String)
   @ Base.MainInclude ./client.jl:489
 [7] top-level scope
   @ REPL[3]:1
in expression starting at /home/ufechner/repos/Tethers.jl/src/Tether_01.jl:21

Any idea what is going on here?

Using collect looks ugly in my eyes. Is there a way to avoid it?

2

If you leave the broadcast out when writing the equations, it works without the collect calls:

eqs = [D(pos) ~ vel
       D(vel) ~ acc
       acc    ~ G_EARTH]
1 Like
3

Thanks a lot! So the first example is fixed, but there are more issues with the latest version of ModelingToolkit.

I have this code:

# Example two: Falling mass, attached to linear spring
using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra

G_EARTH::Vector{Float64} = [0.0, 0.0, -9.81]    # gravitational acceleration     [m/sĀ²]
L0::Float64 = -10.0                             # initial spring length      [m]

# defing the model, Z component upwards
@parameters mass=1.0 c_spring=50.0 damping=0.5 l0=L0
@independent_variables t
@variables   pos(t)[1:3] = [0.0, 0.0,  L0]
@variables   vel(t)[1:3] = [0.0, 0.0,  0.0] 
@variables   acc(t)[1:3] = G_EARTH
@variables unit_vector(t)[1:3]  = [0.0, 0.0, -sign(L0)]
@variables spring_force(t)[1:3] = [0.0, 0.0, 0.0]
@variables norm1(t) = abs(l0) spring_vel(t) = 0.0
D = Differential(t)

eqs = vcat(D(pos)      ~ vel,
           D(vel)      ~ acc,
           norm1        ~ norm(pos),
           unit_vector  ~ -pos/norm1,         # direction from point mass to origin
           spring_vel   ~ -unit_vector ā‹… vel,
           spring_force ~ (c_spring * (norm1 - abs(l0)) + damping * spring_vel) * unit_vector,
           acc          ~ G_EARTH + spring_force/mass)

@named sys = ODESystem(eqs, t)
simple_sys = structural_simplify(sys)

# running the simulation
duration = 10.0
dt = 0.02
tol = 1e-6
tspan = (0.0, duration)
ts    = 0:dt:duration

prob = ODEProblem(simple_sys, nothing, tspan)
@time sol = solve(prob, Rodas5(), dt=dt, abstol=tol, reltol=tol, saveat=ts)

nothing

It works fine with older versions of ModelingToolkit, with the latest version I get this output:

julia> include("src/Tether_02.jl")
ā”Œ Warning: Internal error: Variable (pos(t))[1] was marked as being in 0 ~ sqrt(Symbolics._mapreduce(#396, +, pos(t), Colon(), (:init => false,))) - norm1(t), but was actually zero
ā”” @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/jEIdO/src/structural_transformation/utils.jl:233
ā”Œ Warning: Internal error: Variable (pos(t))[2] was marked as being in 0 ~ sqrt(Symbolics._mapreduce(#396, +, pos(t), Colon(), (:init => false,))) - norm1(t), but was actually zero
ā”” @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/jEIdO/src/structural_transformation/utils.jl:233
ā”Œ Warning: Internal error: Variable (pos(t))[3] was marked as being in 0 ~ sqrt(Symbolics._mapreduce(#396, +, pos(t), Colon(), (:init => false,))) - norm1(t), but was actually zero
ā”” @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/jEIdO/src/structural_transformation/utils.jl:233
ā”Œ Warning: Internal error: Variable (unit_vector(t))[1] was marked as being in 0 ~ -spring_vel(t) + (broadcast(-, unit_vector(t)))[1]*(vel(t))[1] + (broadcast(-, unit_vector(t)))[2]*(vel(t))[2] + (broadcast(-, unit_vector(t)))[3]*(vel(t))[3], but was actually zero
ā”” @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/jEIdO/src/structural_transformation/utils.jl:233
ā”Œ Warning: Internal error: Variable (unit_vector(t))[2] was marked as being in 0 ~ -spring_vel(t) + (broadcast(-, unit_vector(t)))[1]*(vel(t))[1] + (broadcast(-, unit_vector(t)))[2]*(vel(t))[2] + (broadcast(-, unit_vector(t)))[3]*(vel(t))[3], but was actually zero
ā”” @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/jEIdO/src/structural_transformation/utils.jl:233
ā”Œ Warning: Internal error: Variable (unit_vector(t))[3] was marked as being in 0 ~ -spring_vel(t) + (broadcast(-, unit_vector(t)))[1]*(vel(t))[1] + (broadcast(-, unit_vector(t)))[2]*(vel(t))[2] + (broadcast(-, unit_vector(t)))[3]*(vel(t))[3], but was actually zero
ā”” @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/jEIdO/src/structural_transformation/utils.jl:233
  0.280586 seconds (727.74 k allocations: 41.242 MiB, 3.02% gc time, 95.96% compilation time)

Any idea? Is this a bug in ModelingToolkit, or is something wrong with my code?

@ChrisRackauckas Any comment?

4

Open an issue for that.

1 Like
5
6

The issue was closed, but the suggested workaround does not help, it even makes things worse.

Can this issue be re-opened?

7

Ok, the new workaround helps:

@named sys = ODESystem(Symbolics.scalarize.(reduce(vcat, Symbolics.scalarize.(eqs))), t)

But it does not look like the final solution, this should not be needed.