ODE usage: passing parameters used for indexing

New to Julia
,
1

Hi all, I’m trying to solve a 2nd order (mass-spring-damper) ODE by splitting it into two first orders.

In the process of accumulating the force vector in my accel() function, I need to access 4 lists, a points list, a connectivity list, an edge list, and a resting-spring-length list. These aren’t actually parameters of the ODE function, but are necessary for accumulating the force vector. This is, maybe the problem, I’m not sure?

Here is the error I’m getting:

ERROR: LoadError: DimensionMismatch("expected axes([1.0, 0.0, 0.0], 1) = Base.OneTo(1)")

Here is the code for the accel() function:

function accel(a, s, params, t) 

    PL, CL, EL, L, m, f, c = params

    # Add pressure forces 

    for el = 1:size(CL,2)

        ea = (PL[:, CL[2,el]]+s[3*CL[2,el].-[2; 1; 0;]]) - (PL[:, CL[1,el]]+s[3*CL[1,el].-[2; 1; 0;]])

        eb = (PL[:, CL[3,el]]+s[3*CL[3,el].-[2; 1; 0;]]) - (PL[:, CL[1,el]]+s[3*CL[1,el].-[2; 1; 0;]])

        fe = (p/6)*cross(ea,eb)

        ne = [3*CL[1,el].-[2; 1; 0] 3*CL[2,el].-[2;1;0] 3*CL[3,el].-[2; 1; 0]]

        f[ne] = repeat(fe,3,1)    

    end

    # Add spring & damping forces  

    for ele = 1:size(EL,2)

        P = [ PL[:, EL[1,ele]]+s[3*EL[1,ele].-[2; 1; 0;]]   PL[:, EL[2,ele]]+s[3*EL[2,ele].-[2; 1; 0;]] ] # Coords of nodes belonging to this spring

        V = [ s[(3*EL[1,ele].-[2; 1; 0;]).+96]     s[(3*EL[2,ele].-[2; 1; 0;]).+96] ] # Velocities of nodes belonging to this spring

        

        dL = sqrt(sum((P[:,2] - P[:,1]).^2)) # Length of deformed spring (after t = 0)

        

        disp = dL - L[ele] # Displacement from rest

        relV = (V[:,1] - V[:,2]) # Relative velocity

        n_hat = (P[:,2] - P[:,1])/dL # Directional unit vector

        

        spr = L[ele]*disp # Spring force

        damp = c[ele]*relV*n_hat # Damping force

          

        ne = [3*EL[1,ele].-[2; 1; 0] 3*EL[2,ele].-[2; 1; 0]] # Indexing for each node within force vector 

        f[ne[:,1]] .+= -(spr + damp)*n_hat # For node1 of edge, f = k*(norm(pt2-pt1) - restLength)*(pt2-pt1)/norm(pt2-pt1)

        f[ne[:,2]] .+= -(spr + damp)*-n_hat # Same for node2 of edge, except direction flipped, so that f = k*(norm(pt2-pt1) - restLength)*(pt1-pt2)/norm(pt2-pt1)

    end

    a[1:96] = s[1:96]

    a[97:192] .= f ./ m

    

    return a

end

And my call:

params = (PL, CL, EL, L, m, f, c)

prob = ODEProblem(accel, s₀, (0.0,10.0), params)

sol = solve(prob)

Is the error generated because I’m passing in those 4 lists, that don’t have the form of my output vector?

Any insight is appreciated, thank you!

2

Passing lists like that is fine, the parameters variable is just passed through the DifferentialEquations machinery unmodified so you can format it however you like.

That error looks like a broadcast between different size arrays, unfortunately I can’t guess where. Would you mind posting the full stack trace?

You could also try calling your accel method yourself with your initial conditions, that might produce an easier to understand stack trace.

du = similar(s₀)

accel(du, s₀, params, 0.0)
1 Like
3

Ok, yeah you’re right, it’s a broadcast error.

Thanks for the tip on calling accel, much easier to read the stack trace!