ModelingToolkit - how to reset integrator

Specific Domains Modelling & Simulations
, ,
1

I have the following code:

# Falling mass, attached to non-linear spring without compression stiffness,
# initially moving upwards with 4 m/s.
using ModelingToolkit, OrdinaryDiffEq, PyPlot, LinearAlgebra

G_EARTH  = Float64[0.0, 0.0, -9.81]    # gravitational acceleration [m/s²]
L0 = -10.0                             # initial spring length      [m]
V0 = 4                                 # initial velocity           [m/s]

# model, Z component upwards
@parameters mass=1.0 c_spring0=50.0 damping=0.5 l0=L0
@variables t pos(t)[1:3] = [0.0, 0.0,  L0]
@variables   vel(t)[1:3] = [0.0, 0.0,  V0] 
@variables   acc(t)[1:3] = G_EARTH
@variables unit_vector(t)[1:3]  = [0.0, 0.0, -sign(L0)]
@variables c_spring(t) = c_spring0
@variables spring_force(t)[1:3] = [0.0, 0.0, 0.0]
@variables force(t) = 0.0 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,
           c_spring     ~ c_spring0 * (norm1 > abs(l0)),
           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)

duration = 10.0
dt = 0.02
tol = 1e-6
tspan = (0.0, duration)
ts    = 0:dt:duration
u0 = zeros(6)
u0[3] = L0
u0[6] = V0

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

X = sol.t
POS_Z = sol(X, idxs=pos[3])
VEL_Z = sol(X, idxs=vel[3])

lns1 = plot(X, POS_Z, color="green", label="pos_z")
xlabel("time [s]")
ylabel("pos_z [m]")
lns2 = plot(X, L0.+0.005 .* sol[c_spring], color="grey", label="c_spring")
grid(true)
legend() 
twinx()
ylabel("vel_z [m/s]") 
lns3 = plot(X, VEL_Z, color="red", label="vel_z")
lns = vcat(lns1, lns2, lns3)
labs = [l.get_label() for l in lns]
legend(lns, labs) 
PyPlot.show(block=false)
nothing

Output is the following plot:
Nonlinear

This is a mass, attached to a non-linear spring-damper element. I was reading that for better accuracy I should reset the integrator when a discontinuity is happening.

How could I do this?

See also: GitHub - ufechner7/Tethers.jl: Tether models

This is material intended for teaching, so if you see anything I could improve please tell me…

1 Like
2

I added a continues callback:

# Example two: Falling mass, attached to non-linear spring without compression stiffness,
# initially moving upwards with 4 m/s.
using ModelingToolkit, OrdinaryDiffEq, PyPlot, LinearAlgebra

G_EARTH  = Float64[0.0, 0.0, -9.81]    # gravitational acceleration [m/s²]
L0::Float64 = -10.0                             # initial spring length      [m]
V0::Float64 = 4                                 # initial velocity           [m/s]

# model, Z component upwards
@parameters mass=1.0 c_spring0=50.0 damping=0.5 l0=L0
@variables t pos(t)[1:3] = [0.0, 0.0,  L0]
@variables   vel(t)[1:3] = [0.0, 0.0,  V0] 
@variables   acc(t)[1:3] = G_EARTH
@variables unit_vector(t)[1:3]  = [0.0, 0.0, -sign(L0)]
@variables c_spring(t) = c_spring0
@variables spring_force(t)[1:3] = [0.0, 0.0, 0.0]
@variables force(t) = 0.0 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,
           c_spring     ~ c_spring0 * (norm1 > abs(l0)),
           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)

duration = 10.0
dt = 0.02
tol = 1e-6
tspan = (0.0, duration)
ts    = 0:dt:duration
u0 = zeros(6)
u0[3] = L0
u0[6] = V0

function condition(u, t, integrator) # Event when condition(u,t,integrator) == 0
    norm(u[1:3]) - abs(L0)
end
function affect!(integrator)
    println(integrator.t)
end
cb = ContinuousCallback(condition, affect!)

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

X = sol.t
POS_Z = sol(X, idxs=pos[3])
VEL_Z = sol(X, idxs=vel[3])

lns1 = plot(X, POS_Z, color="green", label="pos_z")
xlabel("time [s]")
ylabel("pos_z [m]")
lns2 = plot(X, L0.+0.005 .* sol[c_spring], color="grey", label="c_spring")
grid(true)
legend() 
twinx()
ylabel("vel_z [m/s]") 
lns3 = plot(X, VEL_Z, color="red", label="vel_z")
lns = vcat(lns1, lns2, lns3)
labs = [l.get_label() for l in lns]
legend(lns, labs) 
PyPlot.show(block=false)
nothing

Output:

julia> include("src/Tether_03b.jl")
0.7666377265722647
1.3273695432982244
1.8901011150817926
2.4866925591958946
2.892980433666367
3.5369462103982885
3.8106688933914388
4.526415924233006
4.664257548211656
  0.745579 seconds (266.79 k allocations: 16.595 MiB, 2.25% gc time, 99.49% compilation time: 100% of which was recompilation)

So it kind of works, but:

  • does it actually reset the integrator?
  • is it needed?

It slows down the code a lot, each time I include this file it recompiles which is not happening with the first example…

3

Yes

Depends on the accuracy requirements. In some cases it can make things faster, in other cases it will slow things down.

1 Like
4

Thank you!

Extra question: Is it possible to use static arrays with ModelingToolkit?
Normally an SVector is faster than a normal array…

5

If you make the initial condition a static array then it’ll use static arrays

1 Like