Merge pull request #392 from SebastianM-C/smc/init

Fix initialization for TranslationalPosition and add more tests

Author: ChrisRackauckas

docs/src/connectors/connections.md

@@ -172,7 +172,7 @@ Now, let's consider the position-based approach.  We can build the same model wi
 const TP = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition
 
 systems = @named begin
-    damping = TP.Damper(d = 1, va = 1, vb = 0.0)
+    damping = TP.Damper(d = 1)
     body = TP.Mass(m = 1, v = 1)
     ground = TP.Fixed(s_0 = 0)
 end
@@ -191,7 +191,7 @@ nothing # hide
 As can be seen, we get exactly the same result.  The only difference here is that we are solving an extra equation, which allows us to plot the body position as well.
 
 ```@example connections
-prob = ODEProblem(sys, [], (0, 10.0), [])
+prob = ODEProblem(sys, [], (0, 10.0), [], fully_determined=true)
 sol_p = solve(prob)
 
 p1 = plot(sol_p, idxs = [body.v])
@@ -226,21 +226,21 @@ In this problem, we have a mass, spring, and damper which are connected to a fix
 
 #### Damper
 
-The damper will connect the flange/flange 1 (`flange_a`) to the mass, and flange/flange 2 (`flange_b`) to the fixed point.  For both position- and velocity-based domains, we set the damping constant `d=1` and `va=1` and leave the default for `v_b_0` at 0.  For the position domain, we also need to set the initial positions for `flange_a` and `flange_b`.
+The damper will connect the flange/flange 1 (`flange_a`) to the mass, and flange/flange 2 (`flange_b`) to the fixed point.  For both position- and velocity-based domains, we set the damping constant `d=1` and `va=1` and leave the default for `v_b_0` at 0.
 
 ```@example connections
 @named dv = TV.Damper(d = 1)
-@named dp = TP.Damper(d = 1, va = 1, vb = 0.0, flange_a__s = 3, flange_b__s = 1)
+@named dp = TP.Damper(d = 1)
 nothing # hide
 ```
 
 #### Spring
 
-The spring will connect the flange/flange 1 (`flange_a`) to the mass, and flange/flange 2 (`flange_b`) to the fixed point.  For both position- and velocity-based domains, we set the spring constant `k=1`.  The velocity domain then requires the initial velocity `va` and initial spring stretch `delta_s`.  The position domain instead needs the initial positions for `flange_a` and `flange_b` and the natural spring length `l`.
+The spring will connect the flange/flange 1 (`flange_a`) to the mass, and flange/flange 2 (`flange_b`) to the fixed point.  For both position- and velocity-based domains, we set the spring constant `k=1`.  The velocity domain then requires the initial velocity `va` and initial spring stretch `delta_s`. The position domain instead needs the natural spring length `l`.
 
 ```@example connections
 @named sv = TV.Spring(k = 1)
-@named sp = TP.Spring(k = 1, flange_a__s = 3, flange_b__s = 1, l = 1)
+@named sp = TP.Spring(k = 1, l = 1)
 nothing # hide
 ```
 
@@ -281,7 +281,7 @@ function simplify_and_solve(damping, spring, body, ground; initialization_eqs =
 
     println.(full_equations(sys))
 
-    prob = ODEProblem(sys, [], (0, 10.0), []; initialization_eqs)
+    prob = ODEProblem(sys, [], (0, 10.0), []; initialization_eqs, fully_determined=true)
     sol = solve(prob; abstol = 1e-9, reltol = 1e-9)
 
     return sol

src/Mechanical/TranslationalPosition/components.jl

@@ -16,7 +16,7 @@ Flange fixed in housing at a given position.
         s_0 = 0
     end
     @components begin
-        flange = Flange(; s = s_0)
+        flange = Flange()
     end
     @equations begin
         flange.s ~ s_0
@@ -53,7 +53,7 @@ Sliding mass with inertia
         f(t)
     end
     @components begin
-        flange = Flange(; s = s)
+        flange = Flange()
     end
     @equations begin
         flange.s ~ s
@@ -65,7 +65,7 @@ end
 
 const REL = Val(:relative)
 @component function Spring(::Val{:relative}; name, k, va = 0.0, vb = 0.0,
-        delta_s = 0, flange_a__s = 0, flange_b__s = 0)
+        delta_s = 0)
     pars = @parameters begin
         k = k
     end
@@ -87,56 +87,47 @@ const REL = Val(:relative)
            flange_b.f ~ -f]
 
     return compose(
-        ODESystem(eqs, t, vars, pars; name = name,
-            defaults = [
-                flange_a.s => flange_a__s,
-                flange_b.s => flange_b__s,
-                flange_a.f => +delta_s * k,
-                flange_b.f => -delta_s * k
-            ]),
+        ODESystem(eqs, t, vars, pars; name = name),
         flange_a,
         flange_b)
 end
 
 const ABS = Val(:absolute)
 
 """
-    Spring(; name, k, flange_a__s = 0, flange_b__s = 0, l=0)
+    Spring(; name, k, l=0)
 
 Linear 1D translational spring
 
 # Parameters:
 
   - `k`: [N/m] Spring constant
   - `l`: Unstretched spring length
-  - `flange_a__s`: [m] Initial value of absolute position of flange_a
-  - `flange_b__s`: [m] Initial value of absolute position of flange_b
 
 # Connectors:
 
   - `flange_a: 1-dim. translational flange on one side of spring`
   - `flange_b: 1-dim. translational flange on opposite side of spring` #default function
 """
-function Spring(; name, k, flange_a__s = 0, flange_b__s = 0, l = 0)
-    Spring(ABS; name, k, flange_a__s, flange_b__s, l)
+function Spring(; name, k, l = 0)
+    Spring(ABS; name, k, l)
 end #default function
 
 @component function Spring(::Val{:absolute};
-        name, k, flange_a__s = 0,
-        flange_b__s = 0, l = 0)
+        name, k, l = 0)
     pars = @parameters begin
         k = k
         l = l
     end
     vars = @variables begin
-        f(t) = k * (flange_a__s - flange_b__s - l)
+        f(t)
     end
 
-    @named flange_a = Flange(; s = flange_a__s, f = k * (flange_a__s - flange_b__s - l))
-    @named flange_b = Flange(; s = flange_b__s, f = -k * (flange_a__s - flange_b__s - l))
+    @named flange_a = Flange()
+    @named flange_b = Flange()
 
     eqs = [
-           #    delta_s ~ flange_a.s - flange_b.s
+           #   delta_s ~ flange_a.s - flange_b.s
            f ~ k * (flange_a.s - flange_b.s - l) #delta_s
            flange_a.f ~ +f
            flange_b.f ~ -f]
@@ -166,12 +157,12 @@ Linear 1D translational damper
     @variables begin
         va(t)
         vb(t)
-        f(t) = +(va - vb) * d
+        f(t)
     end
 
     @components begin
-        flange_a = Flange(; s = 0.0, f = (va - vb) * d)
-        flange_b = Flange(; s = 0.0, f = -(va - vb) * d)
+        flange_a = Flange()
+        flange_b = Flange()
     end
 
     @equations begin

test/Mechanical/translational.jl

@@ -32,10 +32,10 @@ end
 
 @testset "Spring, Damper, Mass, Fixed" begin
     @named dv = TV.Damper(d = 1)
-    @named dp = TP.Damper(d = 1, va = 1, vb = 0.0, flange_a.s = 3, flange_b.s = 1)
+    @named dp = TP.Damper(d = 1)
 
     @named sv = TV.Spring(k = 1)
-    @named sp = TP.Spring(k = 1, flange_a__s = 3, flange_b__s = 1, l = 1)
+    @named sp = TP.Spring(k = 1, l = 1)
 
     @named bv = TV.Mass(m = 1)
     @named bp = TP.Mass(m = 1, v = 1, s = 3)
@@ -52,7 +52,8 @@ end
 
         sys = structural_simplify(model)
 
-        prob = ODEProblem(sys, [], (0, 20.0), []; initialization_eqs)
+        prob = ODEProblem(
+            sys, [], (0, 20.0), []; initialization_eqs, fully_determined = true)
         sol = solve(prob; abstol = 1e-9, reltol = 1e-9)
 
         return sol
@@ -71,10 +72,10 @@ end
 
 @testset "driven spring damper mass" begin
     @named dv = TV.Damper(d = 1)
-    @named dp = TP.Damper(d = 1, va = 1.0, vb = 0.0, flange_a.s = 3, flange_b.s = 1)
+    @named dp = TP.Damper(d = 1)
 
     @named sv = TV.Spring(k = 1)
-    @named sp = TP.Spring(k = 1, flange_a__s = 3, flange_b__s = 1, l = 1)
+    @named sp = TP.Spring(k = 1, l = 1)
 
     @named bv = TV.Mass(m = 1)
     @named bp = TP.Mass(m = 1, v = 1, s = 3)
@@ -101,12 +102,13 @@ end
 
     model = System(dv, sv, bv, gv, fv, source)
     sys = structural_simplify(model)
-    prob = ODEProblem(sys, [bv.s => 0], (0, 20.0), [])
+    prob = ODEProblem(
+        sys, [bv.s => 0, sv.delta_s => 1], (0, 20.0), [], fully_determined = true)
     solv = solve(prob, Rodas4())
 
     model = System(dp, sp, bp, gp, fp, source)
     sys = structural_simplify(model)
-    prob = ODEProblem(sys, [], (0, 20.0), [])
+    prob = ODEProblem(sys, [], (0, 20.0), [], fully_determined = true)
     solp = solve(prob, Rodas4())
 
     for sol in (solv, solp)
@@ -186,8 +188,7 @@ end
             function mass_spring(; name)
                 systems = @named begin
                     fixed = TP.Fixed()
-                    spring = TP.Spring(;
-                        k = 10.0, l = 1.0, flange_a__s = 0.0, flange_b__s = 2.0)
+                    spring = TP.Spring(; k = 10.0, l = 1.0)
                     mass = TP.Mass(; m = 100.0, s = 2.0, v = 0.0)
                     pos_sensor = TP.PositionSensor()
                     force_sensor = TP.ForceSensor()
@@ -207,7 +208,7 @@ end
             @named model = mass_spring()
             sys = structural_simplify(model)
 
-            prob = ODEProblem(sys, [], (0.0, 1.0))
+            prob = ODEProblem(sys, [], (0.0, 1.0), fully_determined = true)
             sol = solve(prob, Tsit5())
 
             @test all(sol[sys.spring.flange_a.f] .== sol[sys.force_value.u])
@@ -230,7 +231,7 @@ end
             @named sys = ODESystem(
                 eqs, t, [], []; systems = [force, source, mass, acc, acc_output])
             s = complete(structural_simplify(sys))
-            prob = ODEProblem(s, [mass.s => 0], (0.0, pi))
+            prob = ODEProblem(s, [mass.s => 0], (0.0, pi), fully_determined = true)
             sol = solve(prob, Tsit5())
             @test sol[sys.acc_output.u] ≈ (sol[sys.mass.f] ./ m)
         end