We read every piece of feedback, and take your input very seriously.
To see all available qualifiers, see our documentation.
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
reset_fsal!
Not sure if MTK or OrdinaryDiffEq, but
using ModelingToolkit, OrdinaryDiffEq @component function ContactForce(; name, surface=nothing) systems = @named begin # frame_a = Frame() # frame_b = Frame() # q = RealInput() # force = Force() end vars = @variables begin q(t) = 1 v(t) = 0 f(t) h(t) hd(t) hdd(t) λ(t)=0 end pars = @parameters begin contact::Int = 0 # discrete.time state variable b1::Bool=false # discrete.time state variable b2::Bool=false # discrete.time state variable b3::Bool=false # discrete.time state variable b4::Bool=false # discrete.time state variable a0 = 100 a1 = 20 a2 = 1e6 linFlag::Bool = false end function affect!(integ, u, p, _) h = integ.u[u.h] hd = integ.u[u.hd] λ = integ.u[u.λ] b1_old = integ.p[p.b1] b2_old = integ.p[p.b2] b3_old = integ.p[p.b3] b4_old = integ.p[p.b4] b1 = h <= 0 b2 = λ < 0 b3 = h > 0 b4 = h <= 0 && hd < 0 contact = integ.p[p.contact] integ.p[p.b1] = b1 integ.p[p.b2] = b2 integ.p[p.b3] = b3 integ.p[p.b4] = b4 edge_b1 = b1 != b1_old edge_b2 = b2 != b2_old edge_b3 = b3 != b3_old edge_b4 = b4 != b4_old if contact == 0 && edge_b1 contact = 1 end if contact == 1 && edge_b2 contact = 2 end if contact == 2 && edge_b3 contact = 0 end if contact == 2 && edge_b4 contact = 1 end integ.p[p.contact] = contact end equations = [ 0 ~ ifelse((contact == 1) | (linFlag == true), hdd + a1*hd + a0*h + a2*h^3, λ) f ~ ifelse(linFlag == true, λ, contact*λ) D(q) ~ v 1 * D(v) ~ -1 * 9.81 + f # Temp h ~ q hd ~ D(h) hdd ~ D(hd) ] continuous_events = [h ~ 0, λ ~ 0] => (affect!, [h, hd, λ], [b1,b2,b3,b4,contact], [b1,b2,b3,b4,contact], nothing) ODESystem(equations, t, vars, pars; name, systems, continuous_events) end @named model = ContactForce() model = complete(model) ssys = structural_simplify(model) prob = ODEProblem(ssys, [], (0.0, 5.0)) sol = solve(prob)
julia> sol = solve(prob) ERROR: UndefRefError: access to undefined reference Stacktrace: [1] getproperty @ ./Base.jl:37 [inlined] [2] reset_fsal! @ ~/.julia/packages/OrdinaryDiffEqCore/4A2vD/src/integrators/integrator_utils.jl:485 [inlined] [3] apply_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}) @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/4A2vD/src/integrators/integrator_utils.jl:410 [4] loopheader!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}) @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/4A2vD/src/integrators/integrator_utils.jl:14 [5] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}) @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/4A2vD/src/solve.jl:547 [6] __solve(::ODEProblem{…}, ::CompositeAlgorithm{…}; kwargs::@Kwargs{…}) @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/4A2vD/src/solve.jl:7 [7] __solve @ ~/.julia/packages/OrdinaryDiffEqCore/4A2vD/src/solve.jl:1 [inlined] [8] solve_call(_prob::ODEProblem{…}, args::CompositeAlgorithm{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…}) @ DiffEqBase ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:612 [9] solve_call @ ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:569 [inlined] [10] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::MTKParameters{…}, args::CompositeAlgorithm{…}; kwargs::@Kwargs{…}) @ DiffEqBase ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:1080 [11] solve_up @ ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:1066 [inlined] [12] #solve#51 @ ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:1003 [inlined] [13] solve @ ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:993 [inlined] [14] #__solve#3 @ ~/.julia/packages/OrdinaryDiffEqDefault/2qVWT/src/default_alg.jl:46 [inlined] [15] __solve @ ~/.julia/packages/OrdinaryDiffEqDefault/2qVWT/src/default_alg.jl:45 [inlined] [16] #__solve#72 @ ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:1394 [inlined] [17] __solve @ ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:1386 [inlined] [18] solve_call(::ODEProblem{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{}) @ DiffEqBase ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:612 [19] solve_call @ ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:569 [inlined] [20] #solve_up#53 @ ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:1072 [inlined] [21] solve_up @ ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:1066 [inlined] [22] solve(::ODEProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{}) @ DiffEqBase ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:1003 [23] solve(::ODEProblem{…}) @ DiffEqBase ~/.julia/packages/DiffEqBase/sCsah/src/solve.jl:993 [24] top-level scope @ REPL[150]:1 Some type information was truncated. Use `show(err)` to see complete types.
This appears to be a problem with the default alg, if I specify Rodas4() as the solver it works out.
Rodas4()
OrdinaryDiffEqCore v1.3.0 ModelingToolkit v9.33.1 and master
The text was updated successfully, but these errors were encountered:
will fix
Sorry, something went wrong.
oscardssmith
No branches or pull requests
Not sure if MTK or OrdinaryDiffEq, but
This appears to be a problem with the default alg, if I specify
Rodas4()
as the solver it works out.The text was updated successfully, but these errors were encountered: