Skip to content

Commit

Permalink
Merge pull request #307 from SciML/zero_defaults
Browse files Browse the repository at this point in the history
Reduce the extra initial conditions on Translational components
  • Loading branch information
ChrisRackauckas committed Jul 22, 2024
2 parents 6d9272a + 282ad40 commit 5b05572
Show file tree
Hide file tree
Showing 8 changed files with 79 additions and 76 deletions.
18 changes: 9 additions & 9 deletions src/Hydraulic/IsothermalCompressible/components.jl
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ end

vars = @variables begin
p(t) = p_int
dm(t) = 0
dm(t)
end

systems = @named begin
Expand Down Expand Up @@ -79,7 +79,7 @@ Variable length internal flow model of the fully developed incompressible flow f

@variables begin
x(t) = length_int
ddm(t) = 0
ddm(t)
end

vars = []
Expand Down Expand Up @@ -242,8 +242,8 @@ Reduces the flow from `port_a` to `port_b` by `n`. Useful for modeling parallel
end

vars = @variables begin
dm_a(t) = 0
dm_b(t) = 0
dm_a(t)
dm_b(t)
end

systems = @named begin
Expand Down Expand Up @@ -373,9 +373,9 @@ end

vars = @variables begin
x(t) = x_int
dx(t) = 0
dx(t)
rho(t) = liquid_density(port)
drho(t) = 0
drho(t)
vol(t) = dead_volume + area * x_int
end

Expand Down Expand Up @@ -416,7 +416,7 @@ Fixed fluid volume.

vars = @variables begin
rho(t) = liquid_density(port)
drho(t) = 0
drho(t)
end

# let
Expand Down Expand Up @@ -717,7 +717,7 @@ end

vars = @variables begin
x(t) = x_int
dx(t) = 0
dx(t)
end

eqs = [D(x) ~ dx
Expand Down Expand Up @@ -822,7 +822,7 @@ end

vars = @variables begin
x(t) = x_int
dx(t) = 0
dx(t)
end

total_length = length_a_int + length_b_int
Expand Down
34 changes: 17 additions & 17 deletions src/Mechanical/MultiBody2D/components.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,35 +9,35 @@
end

@variables begin
(A(t) = 0), [state_priority = 10]
(dA(t) = 0), [state_priority = 10]
(ddA(t) = 0), [state_priority = 10]
(A(t)), [state_priority = 10]
(dA(t)), [state_priority = 10]
(ddA(t)), [state_priority = 10]

fx1(t) = 0
fy1(t) = 0
fx1(t)
fy1(t)

fx2(t) = 0
fy2(t) = 0
fx2(t)
fy2(t)

x1(t) = x1_0
dx1(t) = 0
dx1(t)

y1(t) = y1_0
dy1(t) = 0
dy1(t)

x2(t) = l + x1_0
dx2(t) = 0
dx2(t)

y2(t) = 0
dy2(t) = 0
y2(t)
dy2(t)

x_cm(t) = l / 2 + x1_0
dx_cm(t) = 0
ddx_cm(t) = 0
dx_cm(t)
ddx_cm(t)

y_cm(t) = 0
dy_cm(t) = 0
ddy_cm(t) = 0
y_cm(t)
dy_cm(t)
ddy_cm(t)
end

@components begin
Expand Down
6 changes: 3 additions & 3 deletions src/Mechanical/Translational/components.jl
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ Sliding mass with inertia

vars = @variables begin
v(t) = v
f(t) = 0
f(t)
end

eqs = [flange.v ~ v
Expand Down Expand Up @@ -123,7 +123,7 @@ end # default
end
vars = @variables begin
delta_s(t) = delta_s
f(t) = 0
f(t)
end

@named flange_a = MechanicalPort(; v = flange_a__v)
Expand All @@ -148,7 +148,7 @@ const ABS = Val(:absolute)
vars = @variables begin
sa(t) = sa
sb(t) = sb
f(t) = 0
f(t)
end

@named flange_a = MechanicalPort(; v = flange_a__v)
Expand Down
2 changes: 1 addition & 1 deletion src/Mechanical/Translational/sources.jl
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ Linear 1D position input source. Set `solves_force=false` to force input force
a = RealInput()
end

vars = @variables v(t) = 0
vars = @variables v(t)

eqs = [v ~ flange.v
D(v) ~ a.u]
Expand Down
6 changes: 3 additions & 3 deletions src/Mechanical/TranslationalModelica/components.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,8 +48,8 @@ Sliding mass with inertia
m = 0.0, [description = "Mass of sliding mass [kg]"]
end
@variables begin
v(t) = 0.0, [description = "Absolute linear velocity of sliding mass [m/s]"]
a(t) = 0.0, [description = "Absolute linear acceleration of sliding mass [m/s^2]"]
v(t), [description = "Absolute linear velocity of sliding mass [m/s]"]
a(t), [description = "Absolute linear acceleration of sliding mass [m/s^2]"]
end
@extend flange_a, flange_b, s = pr = PartialRigid(; L = 0.0, s)
@equations begin
Expand Down Expand Up @@ -106,7 +106,7 @@ Linear 1D translational damper
d = 0.0, [description = "Damping constant [Ns/m]"]
end
@variables begin
lossPower(t) = 0.0, [description = "Power dissipated by the damper [W]"]
lossPower(t), [description = "Power dissipated by the damper [W]"]
end
@equations begin
f ~ d * v_rel
Expand Down
12 changes: 6 additions & 6 deletions src/Mechanical/TranslationalModelica/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -46,8 +46,8 @@ Partial model for the compliant connection of two translational 1-dim. flanges.
@mtkmodel PartialCompliant begin
@extend (flange_a, flange_b) = pt = PartialTwoFlanges()
@variables begin
s_rel(t) = 0.0, [description = "Relative distance between flanges"]
f(t) = 0.0, [description = "Force between flanges"]
s_rel(t), [description = "Relative distance between flanges", guess = 0.0]
f(t), [description = "Force between flanges", guess = 0.0]
end

@equations begin
Expand All @@ -71,9 +71,9 @@ Partial model for the compliant connection of two translational 1-dim. flanges.
@mtkmodel PartialCompliantWithRelativeStates begin
@extend flange_a, flange_b = pt = PartialTwoFlanges()
@variables begin
s_rel(t) = 0.0, [description = "Relative distance between flanges"]
v_rel(t) = 0.0, [description = "Relative linear velocity))"]
f(t) = 0.0, [description = "Forces between flanges"]
s_rel(t), [description = "Relative distance between flanges"]
v_rel(t), [description = "Relative linear velocity))"]
f(t), [description = "Forces between flanges"]
end

@equations begin
Expand Down Expand Up @@ -149,7 +149,7 @@ end
@mtkmodel PartialRigid begin
@extend flange_a, flange_b = ptf = PartialTwoFlanges()
@variables begin
s(t) = 0.0, [description = "Absolute position of center of component"]
s(t), [description = "Absolute position of center of component"]
end
@parameters begin
L = 0.0, [description = "Length of component, from left flange to right flange"]
Expand Down
2 changes: 1 addition & 1 deletion src/Thermal/utils.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
@connector HeatPort begin
T(t) = 273.15 + 20.0
Q_flow(t) = 0.0, [connect = Flow]
Q_flow(t), [connect = Flow]
end
Base.@doc """
HeatPort(; name, T = 273.15 + 20.0, Q_flow = 0.0)
Expand Down
75 changes: 39 additions & 36 deletions test/Mechanical/translational_modelica.jl
Original file line number Diff line number Diff line change
@@ -1,55 +1,58 @@
using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

using ModelingToolkitStandardLibrary.Blocks
import ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica as TP
using ModelingToolkitStandardLibrary.Blocks: Sine
using ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica: Damper, Spring, Mass,
Fixed, Force

@testset "spring damper mass fixed" begin
@named damper = TP.Damper(; d = 1)
@named spring = TP.Spring(; c = 1, s_rel0 = 1)
@named mass = TP.Mass(; m = 1, v = 1)
@named fixed = TP.Fixed(s0 = 1)

eqs = [connect(spring.flange_a, mass.flange_a, damper.flange_a)
connect(spring.flange_b, damper.flange_b, fixed.flange)]

@named model = ODESystem(eqs, t; systems = [fixed, mass, spring, damper])

sys = structural_simplify(model)

foreach(println, full_equations(sys))
@mtkmodel SpringDamperMassFixed begin
@components begin
damper = Damper(; d = 1)
spring = Spring(; c = 1, s_rel0 = 1)
mass = Mass(; m = 1, v = 1, s = 0)
fixed = Fixed(s0 = 1)
end
@equations begin
connect(spring.flange_a, mass.flange_a, damper.flange_a)
connect(spring.flange_b, damper.flange_b, fixed.flange)
end
end

@mtkbuild sys = SpringDamperMassFixed()

prob = ODEProblem(sys, [], (0, 20.0), [])
sol = solve(prob, ImplicitMidpoint(), dt = 0.01)

@test sol[mass.v][1] == 1.0
@test sol[mass.v][end]0.0 atol=1e-4
@test sol[sys.mass.v][1] == 1.0
@test sol[sys.mass.v][end]0.0 atol=1e-4
end

@testset "driven spring damper mass" begin
@named damper = TP.Damper(; d = 1)
@named spring = TP.Spring(; c = 1, s_rel0 = 1)
@named mass = TP.Mass(; m = 1, v = 1)
@named fixed = TP.Fixed(; s0 = 1)
@named force = TP.Force()
@named source = Sine(frequency = 3, amplitude = 2)

eqs = [connect(force.f, source.output)
connect(force.flange, mass.flange_a)
connect(spring.flange_a, mass.flange_b, damper.flange_a)
connect(spring.flange_b, damper.flange_b, fixed.flange)]

@named model = ODESystem(eqs, t;
systems = [fixed, mass, spring, damper, force, source])

sys = structural_simplify(model)

foreach(println, full_equations(sys))
@mtkmodel DrivenSpringDamperMass begin
@components begin
damper = Damper(; d = 1)
spring = Spring(; c = 1, s_rel0 = 1)
mass = Mass(; m = 1, v = 1, s = 0)
fixed = Fixed(; s0 = 1)
force = Force()
source = Sine(frequency = 3, amplitude = 2)
end

@equations begin
connect(force.f, source.output)
connect(force.flange, mass.flange_a)
connect(spring.flange_a, mass.flange_b, damper.flange_a)
connect(spring.flange_b, damper.flange_b, fixed.flange)
end
end

@mtkbuild sys = DrivenSpringDamperMass()

prob = ODEProblem(sys, [], (0, 20.0), [])
sol = solve(prob, Rodas4())

lb, ub = extrema(sol(15:0.05:20, idxs = mass.v).u)
lb, ub = extrema(sol(15:0.05:20, idxs = sys.mass.v).u)
@test -lbub atol=1e-2
@test -0.11 < lb < -0.1
end

0 comments on commit 5b05572

Please sign in to comment.