# ModelingToolkit - how to reset integrator
**URL:** https://discourse.julialang.org/t/modelingtoolkit-how-to-reset-integrator/107364
**Category:** Modelling & Simulations
**Tags:** question, package, modelingtoolkit
**Created:** [December 9, 2023, 5:07pm UTC](https://discourse.julialang.org/t/modelingtoolkit-how-to-reset-integrator/107364 "2023-12-09T17:07:21Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [December 9, 2023, 5:07pm UTC](https://discourse.julialang.org/t/modelingtoolkit-how-to-reset-integrator/107364/1 "2023-12-09T17:07:21Z")
I have the following code:
```julia
# 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:
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](https://github.com/ufechner7/Tethers.jl)
This is material intended for teaching, so if you see anything I could improve please tell me…
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [December 9, 2023, 6:57pm UTC](https://discourse.julialang.org/t/modelingtoolkit-how-to-reset-integrator/107364/2 "2023-12-09T18:57:41Z")
I added a continues callback:
```julia
# 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
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…
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 9, 2023, 8:20pm UTC](https://discourse.julialang.org/t/modelingtoolkit-how-to-reset-integrator/107364/3 "2023-12-09T20:20:12Z")
> [@ufechner7](#):
>
> - does it actually reset the integrator?
Yes
> [@ufechner7](#):
>
> - is it needed?
Depends on the accuracy requirements. In some cases it can make things faster, in other cases it will slow things down.
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [December 9, 2023, 8:52pm UTC](https://discourse.julialang.org/t/modelingtoolkit-how-to-reset-integrator/107364/4 "2023-12-09T20:52:47Z")
Thank you!
Extra question: Is it possible to use static arrays with ModelingToolkit?
Normally an SVector is faster than a normal array…
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 9, 2023, 9:17pm UTC](https://discourse.julialang.org/t/modelingtoolkit-how-to-reset-integrator/107364/5 "2023-12-09T21:17:32Z")
If you make the initial condition a static array then it’ll use static arrays