diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl index 0a9a2f7f2..a3020cbdb 100644 --- a/src/Blocks/sources.jl +++ b/src/Blocks/sources.jl @@ -442,18 +442,7 @@ end # - SawTooth Generate saw tooth signal # - Trapezoid Generate trapezoidal signal of type Real -function linear_interpolation(x1::T, x2::T, t1::T, t2::T, t) where {T <: Real} - if t1 != t2 - slope = (x2 - x1) / (t2 - t1) - intercept = x1 - slope * t1 - - return slope * t + intercept - else - @assert x1==x2 "x1 ($x1) and x2 ($x2) should be equal if t1 == t2" - - return x2 - end -end +# SampledData Parameter struct ---------------- struct Parameter{T <: Real} data::Vector{T} @@ -531,12 +520,53 @@ function Base.show(io::IO, m::MIME"text/plain", p::Parameter) end end -function get_sampled_data(t, memory::Parameter{T}) where {T} +get_sample_time(memory::Parameter) = memory.ref +Symbolics.@register_symbolic get_sample_time(memory) + +Base.convert(::Type{T}, x::Parameter{T}) where {T <: Real} = x.ref +function Base.convert(::Type{<:Parameter{T}}, x::Number) where {T <: Real} + Parameter{T}(T[], x, true) +end + +# SampledData utilities ---------------- + +function linear_interpolation(x1::Real, x2::Real, t1::Real, t2::Real, t) + if t1 != t2 + slope = (x2 - x1) / (t2 - t1) + intercept = x1 - slope * t1 + + return slope * t + intercept + else + @assert x1==x2 "x1 ($x1) and x2 ($x2) should be equal if t1 == t2" + + return x2 + end +end + +function first_order_backwards_difference(t, memory) + Δt = get_sample_time(memory) + x1 = get_sampled_data(t, memory) + x0 = get_sampled_data(t - Δt, memory) + + return (x1 - x0) / Δt +end + +function first_order_backwards_difference(t, buffer, Δt, circular_buffer) + x1 = get_sampled_data(t, buffer, Δt, circular_buffer) + x0 = get_sampled_data(t - Δt, buffer, Δt, circular_buffer) + + return (x1 - x0) / Δt +end + +function get_sampled_data(t, + buffer::Vector{T}, + dt::T, + circular_buffer = true) where {T <: Real} if t < 0 t = zero(t) end - if isempty(memory.data) + if isempty(buffer) if T <: AbstractFloat return T(NaN) else @@ -544,18 +574,18 @@ function get_sampled_data(t, memory::Parameter{T}) where {T} end end - i1 = floor(Int, t / memory.ref) + 1 #expensive + i1 = floor(Int, t / dt) + 1 #expensive i2 = i1 + 1 - t1 = (i1 - 1) * memory.ref - x1 = @inbounds memory.data[i1] + t1 = (i1 - 1) * dt + x1 = @inbounds buffer[i1] if t == t1 return x1 else - n = length(memory.data) + n = length(buffer) - if memory.circular_buffer + if circular_buffer i1 = (i1 - 1) % n + 1 i2 = (i2 - 1) % n + 1 else @@ -565,122 +595,133 @@ function get_sampled_data(t, memory::Parameter{T}) where {T} end end - t2 = (i2 - 1) * memory.ref - x2 = @inbounds memory.data[i2] + t2 = (i2 - 1) * dt + x2 = @inbounds buffer[i2] return linear_interpolation(x1, x2, t1, t2, t) end end - -get_sample_time(memory::Parameter) = memory.ref -Symbolics.@register_symbolic get_sample_time(memory) - -Symbolics.@register_symbolic get_sampled_data(t, memory) - -function first_order_backwards_difference(t, memory) - Δt = get_sample_time(memory) - x1 = get_sampled_data(t, memory) - x0 = get_sampled_data(t - Δt, memory) - - return (x1 - x0) / Δt +function get_sampled_data(t, buffer) + get_sampled_data(t, buffer.data, buffer.ref, buffer.circular_buffer) end +Symbolics.@register_symbolic get_sampled_data(t, buffer) +Symbolics.@register_symbolic get_sampled_data(t, buffer, dt, circular_buffer) false function Symbolics.derivative(::typeof(get_sampled_data), args::NTuple{2, Any}, ::Val{1}) t = @inbounds args[1] - memory = @inbounds args[2] - first_order_backwards_difference(t, memory) + buffer = @inbounds args[2] + first_order_backwards_difference(t, buffer) +end +function ChainRulesCore.frule((_, ẋ, _), ::typeof(get_sampled_data), t, buffer) + first_order_backwards_difference(t, buffer) * ẋ +end + +function Symbolics.derivative(::typeof(get_sampled_data), args::NTuple{4, Any}, ::Val{1}) + t = @inbounds args[1] + buffer = @inbounds args[2] + sample_time = @inbounds args[3] + circular_buffer = @inbounds args[4] + first_order_backwards_difference(t, buffer, sample_time, circular_buffer) +end +function ChainRulesCore.frule((_, ẋ, _), + ::typeof(get_sampled_data), + t, + buffer, + sample_time, + circular_buffer) + first_order_backwards_difference(t, buffer, sample_time, circular_buffer) * ẋ end -function ChainRulesCore.frule((_, ẋ, _), ::typeof(get_sampled_data), t, memory) - first_order_backwards_difference(t, memory) * ẋ + +# SampledData component ---------------- + +module SampledDataType +@enum Option vector_based struct_based end """ - SampledData(; name, buffer) + SampledData(; name, buffer, sample_time, circular_buffer=true) data input component. # Parameters: - - `buffer`: a `Parameter` type which holds the data and sample time + - `buffer::Vector{Real}`: holds the data sampled at `sample_time` + - `sample_time::Real` + - `circular_buffer::Bool = true`: how to handle `t > length(buffer)*sample_time`. If true data is considered circular, otherwise last data point is held. # Connectors: - `output` """ -@component function SampledData(; name, buffer) +@component function SampledData(::Val{SampledDataType.vector_based}; + name, + buffer, + sample_time, + circular_buffer = true) pars = @parameters begin - buffer = buffer + buffer = buffer #::Vector{Real} + sample_time = sample_time #::Real + circular_buffer = circular_buffer #::Bool end vars = [] systems = @named begin output = RealOutput() end eqs = [ - output.u ~ get_sampled_data(t, buffer), + output.u ~ get_sampled_data(t, buffer, sample_time, circular_buffer), ] return ODESystem(eqs, t, vars, pars; name, systems, - defaults = [output.u => get_sampled_data(0.0, buffer)]) + defaults = [ + output.u => get_sampled_data(0.0, buffer, sample_time, circular_buffer), + ]) end -@deprecate Input SampledData -function SampledData(T::Type, circular_buffer = true; name) - SampledData(T[], zero(T), circular_buffer; name) -end -function SampledData(dt::T, circular_buffer = true) where {T <: Real} - SampledData(T[], dt, circular_buffer; name) -end -function SampledData(data::Vector{T}, dt::T, circular_buffer = true; name) where {T <: Real} - SampledData(; name, buffer = Parameter(data, dt, circular_buffer)) -end +""" + SampledData(; name, buffer) -Base.convert(::Type{T}, x::Parameter{T}) where {T <: Real} = x.ref -function Base.convert(::Type{<:Parameter{T}}, x::Number) where {T <: Real} - Parameter{T}(T[], x, true) -end +data input component. -# 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 +# Parameters: + - `buffer`: a `Parameter` type which holds the data and sample time - if t == zero(t) - empty!(memory.data) +# Connectors: + - `output` +""" +@component function SampledData(::Val{SampledDataType.struct_based}; name, buffer) + pars = @parameters begin + buffer = buffer #::Parameter end - - n = length(memory.data) - i = round(Int, t / Δt) + 1 #expensive - if i == n + 1 - push!(memory.data, DiffEqBase.value(x)) - elseif i <= n - @inbounds memory.data[i] = DiffEqBase.value(x) - else - error("Memory buffer skipped a step: n=$n, i=$i") + vars = [] + systems = @named begin + output = RealOutput() 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 ChainRulesCore.frule((_, _, ṫ, ẋ, _), - ::typeof(set_sampled_data!), - memory, - t, - x, - Δt) - first_order_backwards_difference(t, x, Δt, memory) * ṫ + ẋ + eqs = [ + output.u ~ get_sampled_data(t, buffer), + ] + return ODESystem(eqs, t, vars, pars; name, systems, + defaults = [output.u => get_sampled_data(0.0, buffer)]) end -function first_order_backwards_difference(t, x, Δt, memory) - x1 = set_sampled_data!(memory, t, x, Δt) - x0 = get_sampled_data(t - Δt, memory) +SampledData(x::SampledDataType.Option; kwargs...) = SampledData(Val(x); kwargs...) - return (x1 - x0) / Δt +# struct_based +function SampledData(T::Type, circular_buffer = true; name) + SampledData(SampledDataType.struct_based; + name, + buffer = Parameter(T[], zero(T), circular_buffer)) +end + +# vector_based +function SampledData(sample_time::T, circular_buffer = true; name) where {T <: Real} + SampledData(SampledDataType.vector_based; + name, + buffer = T[], + sample_time, + circular_buffer) +end +function SampledData(buffer::Vector{<:Real}, + sample_time::Real, + circular_buffer = true; + name) + SampledData(SampledDataType.vector_based; name, buffer, sample_time, circular_buffer) +end +function SampledData(; name, buffer, sample_time, circular_buffer) + SampledData(SampledDataType.vector_based; name, buffer, sample_time, circular_buffer) end diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index d5ded6358..0ceab5f58 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -116,7 +116,8 @@ Variable length internal flow model of the fully developed incompressible flow f 0 end - eqs = [0 ~ port_a.dm + port_b.dm] + eqs = [0 ~ port_a.dm + port_b.dm + domain_connect(port_a, port_b)] if variable_length push!(eqs, Δp ~ ifelse(c > 0, shear + inertia, zero(c))) @@ -302,6 +303,7 @@ end end eqs = [0 ~ port_a.dm + port_b.dm + domain_connect(port_a, port_b) dm ~ regRoot(2 * Δp * ρ / c) * x y ~ x] @@ -549,7 +551,7 @@ dm ────► │ │ area p_int, x_int = 0, area, - dead_volume = N == 0 ? area * x_max : 0, + dead_volume = N == 0 ? area * x_int : 0, Χ1 = N == 0 ? 1 : 0, Χ2 = 1) diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index b2df09461..b1df8f21b 100644 --- a/test/Blocks/sources.jl +++ b/test/Blocks/sources.jl @@ -417,27 +417,60 @@ 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([y ~ src.output.u - D(y) ~ dy - D(dy) ~ ddy - connect(src.output, int.input)], - t, - systems = [int, src]) - sys = structural_simplify(iosys) - s = complete(iosys) - prob = ODEProblem(sys, [], (0.0, t_end), [s.src.buffer => Parameter(x, dt)]; - tofloat = false) - prob = remake(prob; p = Parameter.(prob.p)) - - sol = solve(prob, Rodas4(); initializealg = NoInit()) - @test sol.retcode == Success - @test sol[src.output.u][1] == 1.0 #check correct initial condition + @testset "using Parameter type" begin + vars = @variables y(t)=1 dy(t)=0 ddy(t)=0 + @named src = SampledData(Float64) + @named int = Integrator() + @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) + s = complete(iosys) + prob = ODEProblem(sys, + [], + (0.0, t_end), + [s.src.buffer => Parameter(x, dt)]; + tofloat = false) + # prob = remake(prob; p = Parameter.(prob.p)) #<-- no longer needed with ModelingToolkit.jl PR #2231 + + sol = solve(prob, Rodas4(); initializealg = NoInit()) + @test sol.retcode == Success + @test sol[src.output.u][1] == 1.0 #check correct initial condition + + @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 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 + @testset "using Vector Based" begin + vars = @variables y(t)=1 dy(t)=0 ddy(t)=0 + @named src = SampledData(dt) + @named int = Integrator() + @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) + s = complete(iosys) + prob = ODEProblem(sys, + [], + (0.0, t_end), + [s.src.buffer => x, s.src.sample_time => dt]; + tofloat = false) + + sol = solve(prob, Rodas4(); initializealg = NoInit()) + @test sol.retcode == Success + @test sol[src.output.u][1] == 1.0 #check correct initial condition + + @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 end diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl index 8ccda0303..2e2e0c836 100644 --- a/test/Hydraulic/isothermal_compressible.jl +++ b/test/Hydraulic/isothermal_compressible.jl @@ -281,6 +281,12 @@ end sys = structural_simplify(system) defs = ModelingToolkit.defaults(sys) s = complete(system) + dt = 1e-4 + time = 0:dt:0.1 + x = @. 0.9 * (time > 0.015) * (time - 0.015)^2 - 25 * (time > 0.02) * (time - 0.02)^3 + defs[s.input.buffer] = Parameter(x, dt) + #prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 0.1); + # tofloat = false)#, jac = true) prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 0.1); tofloat = false, jac = true) @@ -292,15 +298,9 @@ end @test Symbol(defs[s.valve.port_r.ρ]) == Symbol(s.fluid.ρ) @test Symbol(defs[s.snk.port.ρ]) == Symbol(s.fluid.ρ) - dt = 1e-4 - time = 0:dt:0.1 - x = @. 0.9 * (time > 0.015) * (time - 0.015)^2 - 25 * (time > 0.02) * (time - 0.02)^3 - - defs[s.input.buffer] = Parameter(x, dt) # defs[s.piston.Cd_reverse] = 0.1 - p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false)) - prob = remake(prob; p, tspan = (0, time[end])) + prob = remake(prob; tspan = (0, time[end])) @time sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt, initializealg = NoInit())