# 2D arrays with ModelingToolkit

**URL:** <https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448>\
**Category:** Modelling & Simulations\
**Tags:** question, package, modelingtoolkit\
**Created:** [December 11, 2023, 2:33pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448 "2023-12-11T14:33:53Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 11, 2023, 2:33pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/1 "2023-12-11T14:33:54Z")

</div>

I have the following code, representing a mass attached to a spring-damper element:

```julia
using ModelingToolkit, OrdinaryDiffEq, 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)

```

Now I want to model n+1 masses connected by n spring-damper elements to simulate a tether.

Any ideas how to achieve this? Can I define 2D arrays for position, velocity and acceleration? Can I differentiate a 2D array? Can I loop over an array when creating the symbolic equations?

Background: [GitHub - ufechner7/Tethers.jl: Tether models](https://github.com/ufechner7/Tethers.jl)

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 11, 2023, 3:03pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/2 "2023-12-11T15:03:47Z")

</div>

OK, found a first hint: [Symbolic arrays · Symbolics.jl](https://symbolics.juliasymbolics.org/dev/manual/arrays/)

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [December 11, 2023, 4:31pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/4 "2023-12-11T16:31:37Z")

</div>

I’m working on seemingly closely related topics, so pay attention to [Broadcast inside of `@mtkmodel` macro causes precompilation failures when put into a package](https://discourse.julialang.org/t/broadcast-inside-of-mtkmodel-macro-causes-precompilation-failures-when-put-into-a-package/107277) if you intend to use the new syntax at some point

---

<div class="post-metadata">

**Author:** ![fabianmueller](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabianmueller/32/202905_2.png) [@fabianmueller](https://discourse.julialang.org/u/fabianmueller)\
**Post date:** [December 12, 2023, 11:01am UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/5 "2023-12-12T11:01:21Z")

</div>

Could you not try to create an array of single spring mass components and then incorporate the connections between them? Something like this?

```julia
using ModelingToolkit, OrdinaryDiffEq, 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
function springmass(;name)
    @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)
    return ODESystem(eqs, t; name=name)
end

n=5
sys_array = []
for i in range(1, n)
    push!(sys_array, springmass(name=Symbol("springmass$i")))
end
for sys in sys_array
    # Connect the springs in a chain
end

```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 12, 2023, 11:47am UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/6 "2023-12-12T11:47:55Z")

</div>

I have now a nicely working example:

```Julia
# Tutorial example simulating a 3D mass-spring system with a nonlinear spring (no spring forces
# for l < l_0) and n tether segments. 
using ModelingToolkit, OrdinaryDiffEq, Plots, LinearAlgebra

G_EARTH = Float64[0.0, 0.0, -9.81] # gravitational acceleration [m/s²]
L0::Float64 = 10.0 # initial segment length [m]
V0::Float64 = 2 # initial velocity of lowest mass [m/s]
segments::Int64 = 4 # number of tether segments [-]
α0 = π/8 # initial angle [rad]
POS0 = zeros(3, segments+1)
VEL0 = zeros(3, segments+1)
ACC0 = zeros(3, segments+1)
SEGMENTS0 = zeros(3, segments) 
UNIT_VECTORS0 = zeros(3, segments)
for i in 1:segments+1
    l0 = -(i-1)*L0
    v0 = (i-1)*V0/segments
    POS0[:, i] .= [sin(α0) * l0, 0, cos(α0) * l0]
    VEL0[:, i] .= [sin(α0) * v0, 0, cos(α0) * v0]
end
for i in 2:segments+1
    ACC0[:, i] .= G_EARTH
end
for i in 1:segments
    UNIT_VECTORS0[:, i] .= [0, 0, 1.0]
    SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
end

# model, Z component upwards
@parameters mass=1.0 c_spring0=50.0 damping=0.5 l_seg=L0
@variables t 
@variables pos(t)[1:3, 1:segments+1] = POS0
@variables vel(t)[1:3, 1:segments+1] = VEL0
@variables acc(t)[1:3, 1:segments+1] = ACC0
@variables segment(t)[1:3, 1:segments] = SEGMENTS0
@variables unit_vector(t)[1:3, 1:segments] = UNIT_VECTORS0
@variables norm1(t)[1:segments] = l_seg * ones(segments)
@variables rel_vel(t)[1:3, 1:segments] = zeros(3, segments)
@variables spring_vel(t)[1:segments] = zeros(segments)
@variables c_spring(t)[1:segments] = c_spring0 * ones(segments)
@variables spring_force(t)[1:3, 1:segments] = zeros(3, segments)
@variables total_force(t)[1:3, 1:segments] = zeros(3, segments)
D = Differential(t)

eqs1 = vcat(D.(pos) ~ vel,
            D.(vel) ~ acc)
eqs2 = []
for i in segments:-1:1
    global eqs2
    eqs2 = vcat(eqs2, segment[:, i] ~ pos[:, i+1] - pos[:, i])
    eqs2 = vcat(eqs2, norm1[i] ~ norm(segment[:, i]))
    eqs2 = vcat(eqs2, unit_vector[:, i] ~ -segment[:, i]/norm1[i])
    eqs2 = vcat(eqs2, rel_vel[:, i] ~ vel[:, i+1] - vel[:, i])
    eqs2 = vcat(eqs2, spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i])
    eqs2 = vcat(eqs2, c_spring[i] ~ c_spring0 * (norm1[i] > l_seg))
    eqs2 = vcat(eqs2, spring_force[:, i] ~ (c_spring[i] * (norm1[i] - l_seg) + damping * spring_vel[i]) * unit_vector[:, i])
    if i == segments
        eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i])
    else
        eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i]- spring_force[:, i+1])
    end
    eqs2 = vcat(eqs2, acc[:, i+1] .~ G_EARTH + total_force[:, i] / mass)
end
eqs2 = vcat(eqs2, acc[:, 1] .~ zeros(3))
eqs = vcat(eqs1..., eqs2)
     
@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 = Float64[]
for i in 1:segments+1
    global u0
    u0 = vcat(u0, POS0[:, i], VEL0[:, i])
end

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

function plot2d(sol, reltime, segments)
    index = Int64(round(reltime*50+1))
    x = Float64[]
    z = Float64[]
    for particle in 1:segments+1
        push!(x, (sol(sol.t, idxs=pos[1, particle]))[index])
        push!(z, (sol(sol.t, idxs=pos[3, particle]))[index])
    end
    x_max = maximum(x)
    z_max = maximum(z)
    plot(x,z, xlabel="x [m]", ylabel="z [m]", legend=false)
    annotate!(15, z_max-3.0, "t=$(round(reltime,digits=1)) s")
    plot!(x, z, seriestype = :scatter) 
    ylims!((-segments*10-10, 0.5))
    xlims!((-segments*5, segments*5))
end

dt = 0.04
for time in 0:dt:10
    display(plot2d(sol, time, segments))
    sleep(0.25*dt)
end
nothing

```

Try it out, it looks cool!

Only one question left:  
I use this code to initialize the state vector of the solver:

```julia
u0 = Float64[]
for i in 1:segments+1
    global u0
    u0 = vcat(u0, POS0[:, i], VEL0[:, i])
end

```

This seams to be a bit fragile. It is not clear to me if the order of the states in the state vector is well defined. Is there a better way to initialize the state vector of the solver?

UPDATE: This is better:

```julia
# initial state
u0 = Dict(pos=>POS0, vel=>VEL0)

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [December 12, 2023, 11:52am UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/7 "2023-12-12T11:52:05Z")

</div>

> [@ufechner7](#):
>
> This seams to be a bit fragile. It is not clear to me if the order of the states in the state vector is well defined.

Don’t do this, the order can (and often do) change with any patch release of MTK.

The initial condition is supposed to be supplied as a dict

```julia
Dict(var => value, ...)

```

since symbolic arrays are a bit buggy, you probably have to do

```julia
Dict([
  collect(array_var1) .=> values1;
  collect(array_var2) .=> values2;
])

```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 12, 2023, 11:57am UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/8 "2023-12-12T11:57:02Z")

</div>

For now I don’t want to use components, just build one large model as I did it in the past numerically without modelingtoolkit. See: [https://github.com/ufechner7/KiteModels.jl/blob/main/src/KPS3.jl](https://github.com/ufechner7/KiteModels.jl/blob/main/src/KPS3.jl)

When this kind of 1 to 1 translation works, than I might re-consider to use components… But in the end I probably need a particle system where the springs are not only between two masses, but in a more complex graph for the bridle of the kite…

So treating the kite as a component separate from the tether probably does not make much sense…

What can make sense is to treat the winch, the aerodynamic and the structural model of the kite separate…

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 12, 2023, 12:52pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/9 "2023-12-12T12:52:08Z")

</div>

This actually works:

```julia
# initial state
u0 = Dict(pos=>POS0, vel=>VEL0)

```

But I could not find that in the documentation of ODEProblem…

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 12, 2023, 3:04pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/10 "2023-12-12T15:04:02Z")

</div>

Is `vcat` the best way to add an equation in a loop?

```julia
eqs = vcat(eqs, total_force[:, i] ~ spring_force[:, i])

```

I also tried `push!()`, but it did not work.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [December 12, 2023, 3:06pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/11 "2023-12-12T15:06:56Z")

</div>

> [@ufechner7](#):
>
> I also tried `push!()`, but it did not work.

`push!` only works when adding a single element. To add multiple element, use `append!`.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 12, 2023, 7:15pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/12 "2023-12-12T19:15:46Z")

</div>

Thank you!

One more question: I use events to detect when the length of a tether segment becomes less than the unstretched length, because then the stiffness becomes zero.

For one segment I can write:

```julia
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) # Not needed, just to show that the callback works
end
cb = ContinuousCallback(condition, affect!)

```

But how can I create events if I have an array of segments, and each of them might reach a length below the unstretched length and must generate an event? Can I create an array of conditions? Or an array of callbacks?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 12, 2023, 10:01pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/13 "2023-12-12T22:01:22Z")

</div>

Guess I need to use a [VectorContinuousCallback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#SciMLBase.VectorContinuousCallback) in this case…

The questions that remains is: I just have state vector `u` … How do I know which indices relate to the `pos` and which to the `vel` vector?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [December 13, 2023, 6:37am UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/14 "2023-12-13T06:37:54Z")

</div>

Have you read the MTK docs on events? [Event Handling and Callback Functions · ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/)

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 13, 2023, 7:22am UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/15 "2023-12-13T07:22:11Z")

</div>

No, I was reading only about events in the documentation of DifferentialEquations, see: [Event Handling and Callback Functions · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/)

Strange that there are no cross-links between the two documentations…

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [December 13, 2023, 8:28am UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/16 "2023-12-13T08:28:45Z")

</div>

> [@ufechner7](#):
>
> Strange that there are no cross-links between the two documentations…

The first two paragraphs of the MTK event docs contains two different links to the DiffEq docs. The DiffEq docs do not point to MTK because MTK is a downstream package, of which there are 100s.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 13, 2023, 9:32am UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/17 "2023-12-13T09:32:23Z")

</div>

OK, I tried:

```julia
# Example two: Falling mass, attached to non-linear spring without compression stiffness,
# initially moving upwards with 4 m/s, using a callback to precisely calculate the discontinuities
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; continuous_events = [norm(pos) ~ abs(L0)])
simple_sys = structural_simplify(sys)

duration = 10.0
dt = 0.02
tol = 1e-6
tspan = (0.0, duration)
ts = 0:dt:duration
# initial state
u0 = Dict(pos=>[0,0,L0], vel=>[0,0,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)
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)
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

```

But this does not work, I get the error:

```julia
julia> include("src/Tether_03b.jl")
ERROR: LoadError: UndefVarError: `pos` not defined
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/SymbolicUtils/ssQsQ/src/code.jl:418 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/Symbolics/5BibL/src/build_function.jl:537 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/SymbolicUtils/ssQsQ/src/code.jl:375 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [5] macro expansion
    @ ./none:0 [inlined]
  [6] generated_callfunc(f::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…}, __args::Vararg{…}) where {…}
    @ ModelingToolkit ./none:0 [inlined]
  [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…}, ::Vector{…}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:150
  [8] #346
    @ ~/.julia/packages/ModelingToolkit/dCa81/src/systems/callbacks.jl:403 [inlined]
  [9] determine_event_occurance
    @ ~/.julia/packages/DiffEqBase/LUZAV/src/callbacks.jl:290 [inlined]
 [10] find_callback_time(integrator::OrdinaryDiffEq.ODEIntegrator{…}, callback::ContinuousCallback{…}, counter::Int64)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/callbacks.jl:397
 [11] macro expansion
    @ ~/.julia/packages/DiffEqBase/LUZAV/src/callbacks.jl:131 [inlined]
 [12] find_first_continuous_callback
    @ ~/.julia/packages/DiffEqBase/LUZAV/src/callbacks.jl:126 [inlined]
 [13] find_first_continuous_callback
    @ ~/.julia/packages/DiffEqBase/LUZAV/src/callbacks.jl:124 [inlined]
 [14] handle_callbacks!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/integrators/integrator_utils.jl:324
 [15] _loopfooter!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/integrators/integrator_utils.jl:247
 [16] loopfooter!
    @ ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/integrators/integrator_utils.jl:200 [inlined]
 [17] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/solve.jl:539
 [18] __solve(::ODEProblem{…}, ::Rodas5{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/solve.jl:6
 [19] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/IsxOE/src/solve.jl:1 [inlined]
 [20] solve_call(_prob::ODEProblem{…}, args::Rodas5{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:561
 [21] solve_call
    @ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:527 [inlined]
 [22] #solve_up#42
    @ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:1010 [inlined]
 [23] solve_up
    @ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:996 [inlined]
 [24] #solve#40
    @ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:933 [inlined]
 [25] top-level scope
    @ ~/repos/Tethers.jl/src/Tether_03b.jl:49
 [26] include(fname::String)
    @ Base.MainInclude ./client.jl:489
 [27] top-level scope
    @ REPL[1]:1
in expression starting at /home/ufechner/repos/Tethers.jl/src/Tether_03b.jl:49
Some type information was truncated. Use `show(err)` to see complete types.

```

Any idea?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 13, 2023, 5:30pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/18 "2023-12-13T17:30:48Z")

</div>

Created a bug report: [continuous\_events fails to recognize vectors · Issue #2383 · SciML/ModelingToolkit.jl · GitHub](https://github.com/SciML/ModelingToolkit.jl/issues/2383)

I see that the newest version of `Symbolics.jl` is not used, but this is due to a restriction in `ModelingToolkit.jl`…

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [December 15, 2023, 2:33pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/19 "2023-12-15T14:33:21Z")

</div>

@ufechner7 Did you find an approach to add multiple tether forces to your `acc` equation ? With a for loop or else ?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [December 15, 2023, 7:48pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/20 "2023-12-15T19:48:37Z")

</div>

Yes, the multi-segment tether simulation works fine: [Examples · Tethers.jl](https://ufechner7.github.io/Tethers.jl/dev/examples/#Multi-segment-tether)

Not yet working is the continues callback for this model, but I know how to implement it.

```julia
for i in segments:-1:1
    global eqs2
    eqs2 = vcat(eqs2, segment[:, i] ~ pos[:, i+1] - pos[:, i])
    eqs2 = vcat(eqs2, norm1[i] ~ norm(segment[:, i]))
    eqs2 = vcat(eqs2, unit_vector[:, i] ~ -segment[:, i]/norm1[i])
    eqs2 = vcat(eqs2, rel_vel[:, i] ~ vel[:, i+1] - vel[:, i])
    eqs2 = vcat(eqs2, spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i])
    eqs2 = vcat(eqs2, c_spring[i] ~ c_spring0 * (norm1[i] > l_seg))
    eqs2 = vcat(eqs2, spring_force[:, i] ~ (c_spring[i] * (norm1[i] - l_seg) + damping * spring_vel[i]) * unit_vector[:, i])
    if i == segments
        eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i])
    else
        eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i]- spring_force[:, i+1])
    end
    eqs2 = vcat(eqs2, acc[:, i+1] .~ G_EARTH + total_force[:, i] / mass)
end

```

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [December 15, 2023, 8:12pm UTC](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448/21 "2023-12-15T20:12:27Z")

</div>

This is working indeed when you have a chain of tethers. What I was wondering about was how to connect multiple tethers to the same point, leading to additional terms in the same equation. A bit like a

```julia
for i in tethers
    total_force .+~ spring_force[i]
end

```

I do not know if that could work

[Next page](https://discourse.julialang.org/t/2d-arrays-with-modelingtoolkit/107448.md?page=2)