Merge pull request #186 from SciML/bgc/ae_hack

`Parameter` type optimizations

Author: YingboMa

src/Blocks/sources.jl

@@ -428,11 +428,10 @@ end
 struct Parameter{T <: Real}
     data::Vector{T}
     ref::T
-    n::Int
 end
 
 function Base.isequal(x::Parameter, y::Parameter)
-    b0 = x.n == y.n
+    b0 = length(x.data) == length(y.data)
     if b0
         b1 = all(x.data .== y.data)
         b2 = x.ref == y.ref
@@ -465,7 +464,7 @@ Base.:^(x::Parameter, y::Parameter) = Base.power_by_squaring(x.ref, y.ref)
 Base.isless(x::Parameter, y::Number) = Base.isless(x.ref, y)
 Base.isless(y::Number, x::Parameter) = Base.isless(y, x.ref)
 
-Base.copy(x::Parameter{T}) where {T} = Parameter{T}(copy(x.data), x.ref, x.n)
+Base.copy(x::Parameter{T}) where {T} = Parameter{T}(copy(x.data), x.ref)
 
 function Base.show(io::IO, m::MIME"text/plain", p::Parameter)
     if !isempty(p.data)
@@ -484,9 +483,8 @@ function Parameter(x::T; tofloat = true) where {T <: Real}
         P = T
     end
 
-    return Parameter(P[], x, 0)
+    return Parameter(P[], x)
 end
-Parameter(x::Vector{T}, dt::T) where {T <: Real} = Parameter(x, dt, length(x))
 
 function get_sampled_data(t, memory::Parameter{T}) where {T}
     if t < 0
@@ -505,18 +503,19 @@ function get_sampled_data(t, memory::Parameter{T}) where {T}
     i2 = i1 + 1
 
     t1 = (i1 - 1) * memory.ref
-    x1 = @inbounds getindex(memory.data, i1)
+    x1 = @inbounds memory.data[i1]
 
     if t == t1
         return x1
     else
-        if i2 > memory.n
-            i2 = memory.n
+        n = length(memory.data)
+        if i2 > n
+            i2 = n
             i1 = i2 - 1
         end
 
         t2 = (i2 - 1) * memory.ref
-        x2 = @inbounds getindex(memory.data, i2)
+        x2 = @inbounds memory.data[i2]
         return linear_interpolation(x1, x2, t1, t2, t)
     end
 end
@@ -535,10 +534,13 @@ function first_order_backwards_difference(t, memory)
 end
 
 function Symbolics.derivative(::typeof(get_sampled_data), args::NTuple{2, Any}, ::Val{1})
-    first_order_backwards_difference(args[1], args[2])
+    t = @inbounds args[1]
+    memory = @inbounds args[2]
+    first_order_backwards_difference(t, memory)
 end
 
 SampledData(T::Type; name) = SampledData(T[], zero(T); name)
+SampledData(dt::T) where {T <: Real} = SampledData(T[], dt; name)
 function SampledData(data::Vector{T}, dt::T; name) where {T <: Real}
     SampledData(; name, buffer = Parameter(data, dt))
 end
@@ -571,3 +573,45 @@ end
 @deprecate Input SampledData
 
 Base.convert(::Type{T}, x::Parameter{T}) where {T <: Real} = x.ref
+
+# Beta Code for potential AE Hack ----------------------
+function set_sampled_data!(memory::Parameter{T}, t, x, Δt::Parameter{T}) where {T}
+    if t < 0
+        t = zero(t)
+    end
+
+    if t == zero(t)
+        empty!(memory.data)
+    end
+
+    n = length(memory.data)
+    i = round(Int, t / Δt) + 1 #expensive
+    if i == n + 1
+        push!(memory.data, x)
+    elseif i <= n
+        @inbounds memory.data[i] = x
+    else
+        error("Memory buffer skipped a step: n=$n, i=$i")
+    end
+
+    # memory.ref = Δt
+
+    return x
+end
+Symbolics.@register_symbolic set_sampled_data!(memory, t, x, Δt)
+
+function Symbolics.derivative(::typeof(set_sampled_data!), args::NTuple{4, Any}, ::Val{2})
+    memory = @inbounds args[1]
+    t = @inbounds args[2]
+    x = @inbounds args[3]
+    Δt = @inbounds args[4]
+    first_order_backwards_difference(t, x, Δt, memory)
+end
+Symbolics.derivative(::typeof(set_sampled_data!), args::NTuple{4, Any}, ::Val{3}) = 1 #set_sampled_data returns x, therefore d/dx (x) = 1
+
+function first_order_backwards_difference(t, x, Δt, memory)
+    x1 = set_sampled_data!(memory, t, x, Δt)
+    x0 = get_sampled_data(t - Δt, memory)
+
+    return (x1 - x0) / Δt
+end

src/Hydraulic/IsothermalCompressible/components.jl

@@ -506,7 +506,7 @@ dm ────►               │  │ area
 
     ports = @named begin
         port = HydraulicPort(; p_int)
-        flange = MechanicalPort()
+        flange = MechanicalPort(; f_int = p_int*area)
         damper = ValveBase(; p_a_int = p_int, p_b_int = p_int, area_int = 1, Cd,
             Cd_reverse, minimum_area)
     end
@@ -528,8 +528,7 @@ dm ────►               │  │ area
     Δx = ParentScope(x_max) / N
     x₀ = ParentScope(x_int)
 
-    @named moving_volume = VolumeBase(; p_int, x_int = 0, area, dead_volume = 0, Χ1 = 0,
-        Χ2 = 1)
+    @named moving_volume = VolumeBase(; p_int, x_int = 0, area, dead_volume = 0, Χ1 = 0, Χ2 = 1)
 
     volumes = []
     for i in 1:N
@@ -582,6 +581,8 @@ dm ────►               │  │ area
         defaults = [flange.v => 0])
 end
 
+
+
 @component function SpoolValve(reversible = false; p_a_int, p_b_int, x_int, Cd, d, name)
     pars = @parameters begin
         p_a_int = p_a_int

src/Mechanical/Translational/components.jl

@@ -1,17 +1,17 @@
 """
     Free(; name)
 
-Use to close a system that has un-connected `MechanicalPort`'s.  
+Use to close a system that has un-connected `MechanicalPort`'s where the force should not be zero (i.e. you want to solve for the force to produce the given movement of the port)
 
 # Connectors:
 
   - `flange`: 1-dim. translational flange
 """
 @component function Free(; name)
     @named flange = MechanicalPort()
-    vars = []
+    vars = @variables f(t) = 0
     eqs = [
-        flange.f ~ 0,
+        flange.f ~ f
     ]
     return compose(ODESystem(eqs, t, vars, []; name, defaults = [flange.v => 0]),
         flange)

src/Mechanical/Translational/sources.jl

@@ -36,7 +36,7 @@ Linear 1D position input source
   - `flange`: 1-dim. translational flange
   - `s`: real input 
 """
-@component function Position(; s_0 = 0, name)
+@component function Position(solves_force=true; s_0 = 0, name)
     systems = @named begin
         flange = MechanicalPort()
         s = RealInput()
@@ -45,8 +45,50 @@ Linear 1D position input source
     pars = @parameters s_0 = s_0
     vars = @variables x(t) = s_0
 
-    eqs = [D(x) ~ flange.v
-        s.u ~ x]
+    eqs = [
+        D(x) ~ flange.v
+        s.u ~ x
+        ]
+
+    !solves_force && push!(eqs, 0 ~ flange.f)
 
     ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0, s.u => s_0])
 end
+
+
+@component function Velocity(solves_force=true; name)
+    systems = @named begin
+        flange = MechanicalPort()
+        v = RealInput()
+    end
+
+    pars = []
+    vars = []
+
+    eqs = [
+        v.u ~ flange.v
+        ]
+
+    !solves_force && push!(eqs, 0 ~ flange.f)
+
+    ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0])
+end
+
+
+@component function Acceleration(solves_force=true; s_0 = 0, name)
+    systems = @named begin
+        flange = MechanicalPort()
+        a = RealInput()
+    end
+
+    pars = []
+    vars = @variables v(t) = 0
+
+    eqs = [
+        v ~ flange.v
+        D(v) ~ a.u]
+
+    !solves_force && push!(eqs, 0 ~ flange.f)
+
+    ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0])
+end

src/Mechanical/Translational/utils.jl

@@ -1,10 +1,13 @@
-@connector function MechanicalPort(; name)
-    pars = []
+@connector function MechanicalPort(; name, f_int=0, v_int=0)
+    pars = @parameters begin
+        f_int = f_int
+        v_int = v_int
+    end
     vars = @variables begin
-        v(t)
+        v(t) = v_int
         f(t), [connect = Flow]
     end
-    ODESystem(Equation[], t, vars, pars, name = name, defaults = [f => 0])
+    ODESystem(Equation[], t, vars, pars; name, defaults=[f=>f_int])
 end
 Base.@doc """
     MechanicalPort(;name)

test/Blocks/sources.jl

@@ -6,6 +6,7 @@ using ModelingToolkitStandardLibrary.Blocks: smooth_sin, smooth_cos, smooth_damp
 using OrdinaryDiffEq: ReturnCode.Success
 
 @parameters t
+D = Differential(t)
 
 @testset "Constant" begin
     @named src = Constant(k = 2)
@@ -413,11 +414,14 @@ end
     time = 0:dt:t_end
     x = @. time^2 + 1.0
 
+    vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
+
     @named src = SampledData(Float64)
     @named int = Integrator()
-    @named iosys = ODESystem([
-            connect(src.output, int.input),
-        ],
+    @named iosys = ODESystem([y ~ src.output.u
+            D(y) ~ dy
+            D(dy) ~ ddy
+            connect(src.output, int.input)],
         t,
         systems = [int, src])
     sys = structural_simplify(iosys)
@@ -432,4 +436,6 @@ end
 
     @test sol(time)[src.output.u]≈x atol=1e-3
     @test sol[int.output.u][end]≈1 / 3 * 10^3 + 10.0 atol=1e-3 # closed-form solution to integral
+    @test sol[dy][end]≈2 * time[end] atol=1e-3
+    @test sol[ddy][end]≈2 atol=1e-3
 end

test/Mechanical/translational.jl

@@ -10,11 +10,14 @@ D = Differential(t)
 @testset "Free" begin
     function System(; name)
         systems = @named begin
-            mass = TV.Mass(; m = 100, g = -10)
+            acc = TV.Acceleration(false)
+            a = Constant(k = -10)
+            mass = TV.Mass(; m = 100)
             free = TV.Free()
         end
 
-        eqs = [connect(mass.flange, free.flange)]
+        eqs = [connect(a.output, acc.a)
+               connect(mass.flange, acc.flange, free.flange)]
 
         ODESystem(eqs, t, [], []; name, systems)
     end
@@ -26,6 +29,7 @@ D = Differential(t)
     sol = solve(prob, Rosenbrock23())
 
     @test sol[s.mass.flange.v][end]≈-0.1 * 10 atol=1e-3
+    @test sol[s.free.f][end] ≈ 100 * 10
 end
 
 @testset "spring damper mass fixed" begin