From 885cedaecd4086da4331d31c63753c5c988ac4a0 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Wed, 30 Oct 2024 13:21:23 -0700 Subject: [PATCH 1/4] Add a test for callback (p)reinitialization from Fredrik --- test/symbolic_events.jl | 42 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index b58d5911f4..2bfb34c75e 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1075,3 +1075,45 @@ end prob = ODEProblem(pend, [x => 1], (0.0, 3.0), guesses = [y => x]) @test all(≈(0.0; atol = 1e-9), solve(prob, Rodas5())[[x, y]][end]) end + +@testset "Contact force (or lack thereof)" begin + @component function ContactForce2(; name, surface = nothing) + vars = @variables begin + q(t) = 1 + v(t) = 0 + f(t) + h(t) + hd(t) + hdd(t) + (λ(t) = 0), [irreducible = true] + end + pars = @parameters begin + contact::Int = 0 # discrete.time state variable + a0 = 100 + a1 = 20 + a2 = 1e6 + end + + equations = [0 ~ ifelse((contact == 1), hdd + a1 * hd + a0 * h + a2 * h^3, λ) + f ~ contact * λ + D(q) ~ v + 1 * D(v) ~ -1 * 9.81 + f + h ~ q + hd ~ D(h) + hdd ~ D(hd)] + + function affect!(integ, u, p, _) + @show integ[u.h] + end + continuous_events = ModelingToolkit.SymbolicContinuousCallback( + [h ~ 0], (affect!, [h], [], [], nothing)) + + ODESystem(equations, t, vars, pars; name, continuous_events) + end + @named model = ContactForce2() + model = complete(model) + ssys = structural_simplify(model) + prob = ODEProblem(ssys, [], (0.0, 5.0)) + sol = solve(prob, Rodas4()) + @test sol[model.h][end] < -20.0 # the object has fallen successfully +end From ff2cf4988533205e0e17a4833abbe9011a6c6ab6 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Thu, 31 Oct 2024 14:20:00 -0700 Subject: [PATCH 2/4] Add check that the event fired exactly once --- test/symbolic_events.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 2bfb34c75e..104070e3b9 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1077,6 +1077,7 @@ end end @testset "Contact force (or lack thereof)" begin + triggered = 0 @component function ContactForce2(; name, surface = nothing) vars = @variables begin q(t) = 1 @@ -1103,7 +1104,7 @@ end hdd ~ D(hd)] function affect!(integ, u, p, _) - @show integ[u.h] + triggered += 1 end continuous_events = ModelingToolkit.SymbolicContinuousCallback( [h ~ 0], (affect!, [h], [], [], nothing)) @@ -1116,4 +1117,5 @@ end prob = ODEProblem(ssys, [], (0.0, 5.0)) sol = solve(prob, Rodas4()) @test sol[model.h][end] < -20.0 # the object has fallen successfully + @test triggered == 1 end From 2b367c8bbee5410f0185418e6b5f8f6845e1e328 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 12 Nov 2024 00:56:53 -0800 Subject: [PATCH 3/4] Add additional example of a system that was getting reset by callbacks --- test/symbolic_events.jl | 46 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 104070e3b9..2ba667dff0 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1076,7 +1076,7 @@ end @test all(≈(0.0; atol = 1e-9), solve(prob, Rodas5())[[x, y]][end]) end -@testset "Contact force (or lack thereof)" begin +@testset "Contact force (or lack thereof); #2994" begin triggered = 0 @component function ContactForce2(; name, surface = nothing) vars = @variables begin @@ -1119,3 +1119,47 @@ end @test sol[model.h][end] < -20.0 # the object has fallen successfully @test triggered == 1 end + +@testset "Inductor converter; #2994" begin + @parameters Rd Rsw C1 L1 V1 Rl + @variables I_L1(t) I_V1(t) v1(t) v2(t) [irreducible = true] v3(t) [irreducible = true] + + eqs = [ + I_L1 + I_V1~ 0 + -I_L1 + v2/Rsw + v2/Rd - v3/Rd~ 0 + C1*D(v3) + v3/Rl - v2/Rd + v3/Rd~ 0 + v1~ 10*sin(2*pi*50*t) + -D(I_L1)*L1 + v1 - v2~ 0] + + function D_on(i, u, p,ctx) + i.ps[p.Rd]=1e-3 + end + function D_off(i, u, p,ctx) + i.ps[p.Rd]=1e3 + end + c = ModelingToolkit.SymbolicContinuousCallback( + [v2-v3~0], (D_on, [v2,v3], [Rd], [], nothing);affect_neg =(D_off, [v2,v3], [Rd], [], nothing), + rootfind = SciMLBase.LeftRootFind) + + @mtkbuild pend = ODESystem(eqs, t;continuous_events=c) + + u0 = [ + I_L1=>0.0, + v3=>0.0, + ] + + g=[ + v1=>0.0, + v2=>0.0, + I_V1=>0.0] + p = [ + Rd => 0.001, + Rsw=>1e6, + C1=>0.01, + L1=>0.001, + V1=>10.0, + Rl=>20.0] + prob = ODEProblem(pend, u0, (0.0, 1.5), p, guesses = g) + sol = solve(prob, Rodas5P(),dtmax=1e-4) + @test all(x -> x[1] < 0.5 || x[2] < 10, sol[[t,v3]]) +end From e83f735c19ee415cf9879c4c5003bebf121cba37 Mon Sep 17 00:00:00 2001 From: Benjamin Chung Date: Tue, 12 Nov 2024 00:57:42 -0800 Subject: [PATCH 4/4] Format --- test/symbolic_events.jl | 54 ++++++++++++++++++++--------------------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 2ba667dff0..01d6ba1b03 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1120,46 +1120,46 @@ end @test triggered == 1 end -@testset "Inductor converter; #2994" begin +@testset "Inductor converter; #2994" begin @parameters Rd Rsw C1 L1 V1 Rl @variables I_L1(t) I_V1(t) v1(t) v2(t) [irreducible = true] v3(t) [irreducible = true] - eqs = [ - I_L1 + I_V1~ 0 - -I_L1 + v2/Rsw + v2/Rd - v3/Rd~ 0 - C1*D(v3) + v3/Rl - v2/Rd + v3/Rd~ 0 - v1~ 10*sin(2*pi*50*t) - -D(I_L1)*L1 + v1 - v2~ 0] + eqs = [I_L1 + I_V1 ~ 0 + -I_L1 + v2 / Rsw + v2 / Rd - v3 / Rd ~ 0 + C1 * D(v3) + v3 / Rl - v2 / Rd + v3 / Rd ~ 0 + v1 ~ 10 * sin(2 * pi * 50 * t) + -D(I_L1) * L1 + v1 - v2 ~ 0] - function D_on(i, u, p,ctx) - i.ps[p.Rd]=1e-3 + function D_on(i, u, p, ctx) + i.ps[p.Rd] = 1e-3 end - function D_off(i, u, p,ctx) - i.ps[p.Rd]=1e3 + function D_off(i, u, p, ctx) + i.ps[p.Rd] = 1e3 end c = ModelingToolkit.SymbolicContinuousCallback( - [v2-v3~0], (D_on, [v2,v3], [Rd], [], nothing);affect_neg =(D_off, [v2,v3], [Rd], [], nothing), + [v2 - v3 ~ 0], (D_on, [v2, v3], [Rd], [], nothing); + affect_neg = (D_off, [v2, v3], [Rd], [], nothing), rootfind = SciMLBase.LeftRootFind) - @mtkbuild pend = ODESystem(eqs, t;continuous_events=c) + @mtkbuild pend = ODESystem(eqs, t; continuous_events = c) u0 = [ - I_L1=>0.0, - v3=>0.0, + I_L1 => 0.0, + v3 => 0.0 ] - g=[ - v1=>0.0, - v2=>0.0, - I_V1=>0.0] - p = [ + g = [ + v1 => 0.0, + v2 => 0.0, + I_V1 => 0.0] + p = [ Rd => 0.001, - Rsw=>1e6, - C1=>0.01, - L1=>0.001, - V1=>10.0, - Rl=>20.0] + Rsw => 1e6, + C1 => 0.01, + L1 => 0.001, + V1 => 10.0, + Rl => 20.0] prob = ODEProblem(pend, u0, (0.0, 1.5), p, guesses = g) - sol = solve(prob, Rodas5P(),dtmax=1e-4) - @test all(x -> x[1] < 0.5 || x[2] < 10, sol[[t,v3]]) + sol = solve(prob, Rodas5P(), dtmax = 1e-4) + @test all(x -> x[1] < 0.5 || x[2] < 10, sol[[t, v3]]) end