diff --git a/Project.toml b/Project.toml
index 8896d1590..d9718057e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,9 +17,10 @@ julia = "1.6"
 
 [extras]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["OrdinaryDiffEq", "SafeTestsets", "Test", "ControlSystemsBase"]
+test = ["OrdinaryDiffEq", "SafeTestsets", "Test", "ControlSystemsBase", "DataInterpolations"]
diff --git a/docs/pages.jl b/docs/pages.jl
index 24ded55eb..05b09ec83 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -5,6 +5,7 @@ pages = [
         "Custom Components" => "tutorials/custom_component.md",
         "Thermal Conduction Model" => "tutorials/thermal_model.md",
         "DC Motor with Speed Controller" => "tutorials/dc_motor_pi.md",
+        "Input Component" => "tutorials/input_component.md",
     ],
     "About Acausal Connections" => "connectors/connections.md",
     "API" => [
diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md
new file mode 100644
index 000000000..8f40a9575
--- /dev/null
+++ b/docs/src/tutorials/input_component.md
@@ -0,0 +1,147 @@
+# Running Models with Discrete Data
+
+There are 2 ways to include data as part of a model.
+
+ 1. using `ModelingToolkitStandardLibrary.Blocks.TimeVaryingFunction`
+ 2. using a custom component with external data or sensor
+ 3. using `ModelingToolkitStandardLibrary.Blocks.Input`
+
+This tutorial demonstrate each case and explain the pros and cons of each.
+
+## TimeVaryingFunction Component
+
+This component is easy to use and is performative.  However the data is locked to the `ODESystem` and can only be changed by building a new `ODESystem`.  Therefore, running a batch of data would not be efficient.  Below is an example of how to use the `TimeVaryingFunction` with `DataInterpolations` to build the function from discrete data.
+
+```julia
+using ModelingToolkit
+using ModelingToolkitStandardLibrary.Blocks
+using DataInterpolations
+using OrdinaryDiffEq
+
+@parameters t
+D = Differential(t)
+
+function System(f; name)
+    src = TimeVaryingFunction(f)
+
+    vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
+    pars = @parameters m=10 k=1000 d=1
+
+    eqs = [f ~ src.output.u
+           ddx * 10 ~ k * x + d * dx + f
+           D(x) ~ dx
+           D(dx) ~ ddx]
+
+    ODESystem(eqs, t, vars, pars; systems = [src], name)
+end
+
+dt = 4e-4
+time = 0:dt:0.1
+data = sin.(2 * pi * time * 100) # example data
+
+f = LinearInterpolation(data, time)
+
+@named system = System(f)
+sys = structural_simplify(system)
+prob = ODEProblem(sys, [], (0, time[end]))
+sol = solve(prob, ImplicitEuler())
+```
+
+If we want to run a new data set, this requires building a new `LinearInterpolation` and `ODESystem` followed by running `structural_simplify`, all of which takes time.  Therefore, to run serveral pieces of data it's better to re-use an `ODESystem`.  The next couple methods will demonstrate this.
+
+## Custom Component with External Data
+
+The below code shows how to include data using a `Ref` and registered `input` function.  This example uses a very basic function which requires non-adaptive solving and sampled data.  As can be seen, the data can easily be set and changed before solving.
+
+```julia
+const rdata = Ref{Vector{Float64}}()
+
+# Data Sets
+data1 = sin.(2 * pi * time * 100)
+data2 = cos.(2 * pi * time * 50)
+
+function input(t)
+    i = floor(Int, t / dt) + 1
+    x = rdata[][i]
+
+    return x
+end
+
+Symbolics.@register_symbolic input(t)
+
+function System(; name)
+    vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
+    pars = @parameters m=10 k=1000 d=1
+
+    eqs = [f ~ input(t)
+           ddx * 10 ~ k * x + d * dx + f
+           D(x) ~ dx
+           D(dx) ~ ddx]
+
+    ODESystem(eqs, t, vars, pars; name)
+end
+
+@named system = System()
+sys = structural_simplify(system)
+prob = ODEProblem(sys, [], (0, time[end]))
+
+rdata[] = data1
+sol1 = solve(prob, ImplicitEuler(); dt, adaptive = false);
+ddx1 = sol1[sys.ddx]
+
+rdata[] = data2
+sol2 = solve(prob, ImplicitEuler(); dt, adaptive = false)
+ddx2 = sol2[sys.ddx]
+```
+
+The drawback of this method is that the solution observables can be linked to the data `Ref`, which means that if the data changes then the observables are no longer valid.  In this case `ddx` is an observable that is derived directly from the data.  Therefore, `sol1[sys.ddx]` is no longer correct after the data is changed for `sol2`.  Additional code could be added to resolve this issue, for example by using a `Ref{Dict}` that could link a parameter of the model to the data source.  This would also be necessary for parallel processing.
+
+```julia
+# the following test will fail
+@test all(ddx1 .== sol1[sys.ddx]) #returns false
+```
+
+## Input Component
+
+To resolve the issues presented above, the `Input` component can be used which allows for a resusable `ODESystem` and self contained data which ensures a solution which remains valid for it's lifetime.  Now it's possible to also parallize the solving.
+
+```julia
+function System(; name)
+    src = Input(Float64)
+
+    vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0
+    pars = @parameters m=10 k=1000 d=1
+
+    eqs = [f ~ src.output.u
+           ddx * 10 ~ k * x + d * dx + f
+           D(x) ~ dx
+           D(dx) ~ ddx]
+
+    ODESystem(eqs, t, vars, pars; systems = [src], name)
+end
+
+@named system = System()
+sys = structural_simplify(system)
+s = complete(system)
+prob = ODEProblem(sys, [], (0, time[end]); tofloat = false)
+defs = ModelingToolkit.defaults(sys)
+
+function get_prob(data)
+    defs[s.src.buffer] = Parameter(data, dt)
+    # ensure p is a uniform type of Vector{Parameter{Float64}} (converting from Vector{Any})
+    p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false))
+    remake(prob; p)
+end
+
+prob1 = get_prob(data1)
+prob2 = get_prob(data2)
+
+sol1 = Ref{ODESolution}()
+sol2 = Ref{ODESolution}()
+@sync begin
+    @async sol1[] = solve(prob1, ImplicitEuler())
+    @async sol2[] = solve(prob2, ImplicitEuler())
+end
+```
+
+Note, in the above example, we can build the system with an empty `Input` component, only setting the expected data type: `@named src = Input(Float64)`.  It's also possible to initialize the component with real data: `@named src = Input(data, dt)`.  Additionally note that before running an `ODEProblem` using the `Input` component that the parameter vector should be a uniform type so that Julia is not slowed down by type instability.  Because the `Input` component contains the `buffer` parameter of type `Parameter`, we must generate the problem using `tofloat=false`.  This will initially give a parameter vector of type `Vector{Any}` with a mix of numbers and `Parameter` type.  We can convert the vector to all `Parameter` type by running `p = Parameter.(p)`.  This will wrap all the single values in a `Parameter` type which will be mathmatically equivalent to a `Number`.
diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl
index 7997034b7..dab81e5f6 100644
--- a/src/Blocks/Blocks.jl
+++ b/src/Blocks/Blocks.jl
@@ -17,7 +17,7 @@ export Log, Log10
 include("math.jl")
 
 export Constant, TimeVaryingFunction, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine,
-       Square, Triangular
+       Square, Triangular, Parameter, Input
 include("sources.jl")
 
 export Limiter, DeadZone, SlewRateLimiter
diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl
index 1cefb5ec1..8ec900db1 100644
--- a/src/Blocks/sources.jl
+++ b/src/Blocks/sources.jl
@@ -411,3 +411,160 @@ end
 # - Pulse   Generate pulse signal of type Real
 # - 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
+
+struct Parameter{T <: Real}
+    data::Vector{T}
+    ref::T
+    n::Int
+end
+
+function Base.isequal(x::Parameter, y::Parameter)
+    b0 = x.n == y.n
+    if b0
+        b1 = all(x.data .== y.data)
+        b2 = x.ref == y.ref
+        return b1 & b2
+    else
+        return false
+    end
+end
+
+Base.:*(x::Number, y::Parameter) = x * y.ref
+Base.:*(y::Parameter, x::Number) = Base.:*(x, y)
+Base.:*(x::Parameter, y::Parameter) = x.ref * y.ref
+
+Base.:/(x::Number, y::Parameter) = x / y.ref
+Base.:/(y::Parameter, x::Number) = y.ref / x
+Base.:/(x::Parameter, y::Parameter) = x.ref / y.ref
+
+Base.:+(x::Number, y::Parameter) = x + y.ref
+Base.:+(y::Parameter, x::Number) = Base.:+(x, y)
+Base.:+(x::Parameter, y::Parameter) = x.ref + y.ref
+
+Base.:-(x::Number, y::Parameter) = x - y.ref
+Base.:-(y::Parameter, x::Number) = y.ref - x
+Base.:-(x::Parameter, y::Parameter) = x.ref - y.ref
+
+Base.:^(x::Number, y::Parameter) = Base.power_by_squaring(x, y.ref)
+Base.:^(y::Parameter, x::Number) = Base.power_by_squaring(y.ref, x)
+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)
+
+function Base.show(io::IO, m::MIME"text/plain", p::Parameter)
+    if !isempty(p.data)
+        print(io, p.data)
+    else
+        print(io, p.ref)
+    end
+end
+
+Parameter(x::Parameter) = x
+function Parameter(x::T; tofloat = true) where {T <: Real}
+    if tofloat
+        x = float(x)
+        P = typeof(x)
+    else
+        P = T
+    end
+
+    return Parameter(P[], x, 0)
+end
+Parameter(x::Vector{T}, dt::T) where {T <: Real} = Parameter(x, dt, length(x))
+
+function input(t, memory::Parameter{T}) where {T}
+    if t < 0
+        t = zero(t)
+    end
+
+    if isempty(memory.data)
+        if T isa Float16
+            return NaN16
+        elseif T isa Float32
+            return NaN32
+        elseif T isa Float64
+            return NaN64
+        else
+            return zero(T)
+        end
+    end
+
+    i1 = floor(Int, t / memory.ref) + 1 #expensive
+    i2 = i1 + 1
+
+    t1 = (i1 - 1) * memory.ref
+    x1 = @inbounds getindex(memory.data, i1)
+
+    if t == t1
+        return x1
+    else
+        if i2 > memory.n
+            i2 = memory.n
+            i1 = i2 - 1
+        end
+
+        t2 = (i2 - 1) * memory.ref
+        x2 = @inbounds getindex(memory.data, 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 input(t, memory)
+
+function first_order_backwards_difference(t, memory)
+    Δt = get_sample_time(memory)
+    x1 = input(t, memory)
+    x0 = input(t - Δt, memory)
+
+    return (x1 - x0) / Δt
+end
+
+function Symbolics.derivative(::typeof(input), args::NTuple{2, Any}, ::Val{1})
+    first_order_backwards_difference(args[1], args[2])
+end
+
+Input(T::Type; name) = Input(T[], zero(T); name)
+function Input(data::Vector{T}, dt::T; name) where {T <: Real}
+    Input(; name, buffer = Parameter(data, dt))
+end
+
+"""
+    Input(; name, buffer)
+
+data input component.  
+
+# Parameters:
+  - `buffer`: a `Parameter` type which holds the data and sample time
+
+# Connectors:
+  - `output`
+"""
+@component function Input(; name, buffer)
+    pars = @parameters begin buffer = buffer end
+    vars = []
+    systems = @named begin output = RealOutput() end
+    eqs = [
+        output.u ~ input(t, buffer),
+    ]
+    return ODESystem(eqs, t, vars, pars; name, systems,
+                     defaults = [output.u => input(0.0, buffer)]) #TODO: get initial value from buffer
+end
diff --git a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl
index f2f7eb2e8..8ac4ffd0e 100644
--- a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl
+++ b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl
@@ -3,9 +3,12 @@ Library to model iso-thermal compressible liquid fluid flow
 """
 module IsothermalCompressible
 
-using ModelingToolkit, Symbolics, IfElse
+using ModelingToolkit, Symbolics
+
 using ...Blocks: RealInput, RealOutput
-using ...Mechanical.Translational: MechanicalPort
+using ...Mechanical.Translational: MechanicalPort, Mass
+
+using IfElse: ifelse
 
 @parameters t
 D = Differential(t)
@@ -13,7 +16,7 @@ D = Differential(t)
 export HydraulicPort, HydraulicFluid
 include("utils.jl")
 
-export Source, InputSource, Cap, Pipe, FixedVolume, DynamicVolume
+export Source, InputSource, Cap, Tube, FixedVolume, DynamicVolume
 include("components.jl")
 
 end
diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl
index 98c90a370..da787dae1 100644
--- a/src/Hydraulic/IsothermalCompressible/components.jl
+++ b/src/Hydraulic/IsothermalCompressible/components.jl
@@ -75,6 +75,22 @@ Caps a hydrualic port to prevent mass flow in or out.
     ODESystem(eqs, t, vars, pars; name, systems)
 end
 
+@component function Open(; p_int, name)
+    pars = @parameters p_int = p_int
+
+    vars = @variables begin
+        p(t) = p_int
+        dm(t) = 0
+    end
+
+    systems = @named begin port = HydraulicPort(; p_int = p_int) end
+
+    eqs = [port.p ~ p
+           port.dm ~ dm]
+
+    ODESystem(eqs, t, vars, pars; name, systems)
+end
+
 """
     FixedVolume(; vol, p_int, name)
 
@@ -96,7 +112,7 @@ Fixed fluid volume.
     systems = @named begin port = HydraulicPort(; p_int) end
 
     vars = @variables begin
-        rho(t) = density(port, p_int)
+        rho(t) = liquid_density(port)
         drho(t) = 0
     end
 
@@ -104,29 +120,29 @@ Fixed fluid volume.
     dm = port.dm
 
     eqs = [D(rho) ~ drho
-           rho ~ density(port, port.p)
+           rho ~ full_density(port)
            dm ~ drho * vol]
 
     ODESystem(eqs, t, vars, pars; name, systems)
 end
 
 """
-    PipeBase(; p_int, area, length, perimeter=2*sqrt(area*pi), shape_factor=64, name)
+    TubeBase(; p_int, area, length, perimeter=2*sqrt(area*pi), shape_factor=64, name)
 
-Pipe segement which models purely the fully developed flow friction, ignoring any compressibility.
+Internal flow model of the fully developed flow friction, ignoring any compressibility.
 
 # Parameters:
-- `p_int`: [Pa] initial pressure (set by `p_int` argument)
-- `area`: [m^2] tube cross sectional area (set by `area` argument)
-- `length`: [m] length of the pipe (set by `length` argument)
-- `perimeter`: [m] perimeter of the pipe cross section (set by optional `perimeter` argument, needed only for non-circular pipes)
+- `p_int`: [Pa] initial pressure
+- `area`: [m^2] tube cross sectional area
+- `length`: [m] length of the pipe
+- `perimeter`: [m] perimeter of the pipe cross section (needed only for non-circular pipes)
 - `Φ`: shape factor, see `friction_factor` function (set by optional `shape_factor` argument, needed only for non-circular pipes).  
 
 # Connectors:
 - `port_a`: hydraulic port
 - `port_b`: hydraulic port
 """
-@component function PipeBase(; p_int, area, length, perimeter = 2 * sqrt(area * pi),
+@component function TubeBase(; p_int, area, length, perimeter = 2 * sqrt(area * pi),
                              shape_factor = 64, name)
     pars = @parameters begin
         p_int = p_int
@@ -149,7 +165,7 @@ Pipe segement which models purely the fully developed flow friction, ignoring an
 
     d_h = 4 * area / perimeter
 
-    ρ = (density(port_a, port_a.p) + density(port_b, port_b.p)) / 2
+    ρ = (full_density(port_a) + full_density(port_b)) / 2
     μ = viscosity(port_a)
 
     f = friction_factor(dm, area, d_h, ρ, μ, Φ)
@@ -162,30 +178,34 @@ Pipe segement which models purely the fully developed flow friction, ignoring an
 end
 
 """
-    Pipe(N; p_int, area, length, perimeter=2*sqrt(area*pi), shape_factor=64, name)
+    Tube(N; p_int, area, length, effective_length=length, perimeter = 2 * sqrt(area * pi), shape_factor = 64, name)
 
-Pipe modeled with `N` segements which models the fully developed flow friction and compressibility.
+Tube modeled with `N` segements which models the fully developed flow friction and compressibility.
 
 # Parameters:
-- `p_int`: [Pa] initial pressure (set by `p_int` argument)
-- `area`: [m^2] tube cross sectional area (set by `area` argument)
-- `length`: [m] length of the pipe (set by `length` argument)
-- `perimeter`: [m] perimeter of the pipe cross section (set by optional `perimeter` argument, needed only for non-circular pipes)
-- `Φ`: shape factor, see `friction_factor` function (set by optional `shape_factor` argument, needed only for non-circular pipes).  
+- `p_int`: [Pa] initial pressure 
+- `area`: [m^2] tube cross sectional area 
+- `length`: [m] real length of the tube 
+- `effective_length`: [m] tube length to account for flow development and other restrictions, used to set the `TubeBase` length which calculates the flow resistance
+- `perimeter`: [m] perimeter of the tube cross section (set by optional `perimeter` argument, needed only for non-circular tubes)
+- `Φ`: shape factor, see `friction_factor` function (set by optional `shape_factor` argument, needed only for non-circular tubes)
 
 # Connectors:
 - `port_a`: hydraulic port
 - `port_b`: hydraulic port
 """
-@component function Pipe(N; p_int, area, length, perimeter = 2 * sqrt(area * pi),
+@component function Tube(N; p_int, area, length, effective_length = length,
+                         perimeter = 2 * sqrt(area * pi),
                          shape_factor = 64, name)
     @assert(N>1,
-            "the pipe component must be defined with more than 1 segment (i.e. N>1), found N=$N")
+            "the Tube component must be defined with more than 1 segment (i.e. N>1), found N=$N")
 
+    #TODO: How to set an assert effective_length >= length ??
     pars = @parameters begin
         p_int = p_int
         area = area
         length = length
+        effective_length = effective_length
         perimeter = perimeter
         Φ = shape_factor
     end
@@ -199,9 +219,9 @@ Pipe modeled with `N` segements which models the fully developed flow friction a
 
     pipe_bases = []
     for i in 1:(N - 1)
-        x = PipeBase(; name = Symbol("p$i"), shape_factor = ParentScope(Φ),
+        x = TubeBase(; name = Symbol("p$i"), shape_factor = ParentScope(Φ),
                      p_int = ParentScope(p_int), area = ParentScope(area),
-                     length = ParentScope(length) / (N - 1),
+                     length = ParentScope(effective_length) / (N - 1),
                      perimeter = ParentScope(perimeter))
         push!(pipe_bases, x)
     end
@@ -225,6 +245,49 @@ Pipe modeled with `N` segements which models the fully developed flow friction a
     ODESystem(eqs, t, vars, pars; name, systems = [ports; pipe_bases; volumes])
 end
 
+"""
+    FlowDivider(;p_int, n, name)
+
+Reduces the flow from `port_a` to `port_b` by `n`.  Useful for modeling parallel tubes efficiently by placing a `FlowDivider` on each end of a tube.
+
+# Parameters:
+- `p_int`: [Pa] initial pressure 
+- `n`: divide flow from `port_a` to `port_b` by `n`
+
+# Connectors:
+- `port_a`: full flow hydraulic port
+- `port_b`: part flow hydraulic port
+"""
+@component function FlowDivider(; p_int, n, name)
+
+    #TODO: assert n >= 1
+
+    pars = @parameters begin
+        n = n
+        p_int = p_int
+    end
+
+    vars = @variables begin
+        dm_a(t) = 0
+        dm_b(t) = 0
+    end
+
+    systems = @named begin
+        port_a = HydraulicPort(; p_int)
+        port_b = HydraulicPort(; p_int)
+        open = Open(; p_int)
+    end
+
+    eqs = [connect(port_a, port_b, open.port)
+           dm_a ~ port_a.dm
+           dm_b ~ dm_a / n
+           open.dm ~ dm_a - dm_b # extra flow dumps into an open port
+           # port_b.dm ~ dm_b # divided flow goes to port_b
+           ]
+
+    ODESystem(eqs, t, vars, pars; name, systems)
+end
+
 """
     DynamicVolume(; p_int, x_int=0, area, dead_volume=0, direction=+1, minimum_volume=0, name)
 
@@ -240,6 +303,13 @@ dm ────►  dead volume  │  │ area
                        │
                        └─► x (= flange.v * direction)
 ```
+# States:
+- `x(t)`: [m] moving wall position
+- `dx(t)`: [m/s] moving wall velocity
+- `rho(t)`: [kg/m^3] density
+- `drho(t)`: [kg/s-m^3] density derivative
+- `vol(t)`: [m^3] volume
+- `p(t)`: [Pa] dynamic pressure
 
 # Parameters:
 - `p_int`: [Pa] initial pressure
@@ -271,20 +341,25 @@ dm ────►  dead volume  │  │ area
     vars = @variables begin
         x(t) = x_int
         dx(t) = 0
-        rho(t) = density(port, p_int)
+        rho(t) = liquid_density(port)
         drho(t) = 0
         vol(t) = dead_volume + area * x_int
-        p(t) = p_int
+        p(t) = p_int # p represents the dynamic pressure
     end
 
-    eqs = [0 ~ IfElse.ifelse(vol >= minimum_volume, p - port.p, port.dm)
+    # let
+    dm = port.dm
+    u = dm / (rho * area)
+
+    eqs = [0 ~ ifelse(vol >= minimum_volume, (p) - (port.p - (1 / 2) * rho * u^2),
+                      dm)
            vol ~ dead_volume + area * x
            D(x) ~ dx
            D(rho) ~ drho
            dx ~ flange.v * direction
-           rho ~ density(port, p)
-           port.dm ~ drho * vol + rho * area * dx
-           flange.f ~ -p * area * direction] #TODO: update to dynamic pressure
+           rho ~ full_density(port)
+           dm ~ drho * vol + rho * area * dx
+           flange.f ~ -p * area * direction]
 
     ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0])
 end
@@ -305,7 +380,7 @@ Valve with area input and discharge coefficient `Cd` defined by https://en.wikip
 - `port_b`: hydraulic port
 - `input`: real input setting the valve `area`.  Note: absolute value taken
 """
-@component function Valve(; p_a_int, p_b_int, area_int, Cd, name)
+@component function Valve(reversible = false; p_a_int, p_b_int, area_int, Cd, name)
     pars = @parameters begin
         p_a_int = p_a_int
         p_b_int = p_b_int
@@ -322,13 +397,136 @@ Valve with area input and discharge coefficient `Cd` defined by https://en.wikip
     vars = []
 
     # let
-    ρ = (density(port_a, port_a.p) + density(port_b, port_b.p)) / 2
+    ρ = (full_density(port_a) + full_density(port_b)) / 2
     Δp = port_a.p - port_b.p
     dm = port_a.dm
-    area = abs(input.u)
+    area = if reversible
+        input.u
+    else
+        ifelse(input.u > 0, input.u, 0)
+    end
 
     eqs = [sign(Δp) * dm ~ sqrt(2 * abs(Δp) * ρ / Cd) * area
            0 ~ port_a.dm + port_b.dm]
 
     ODESystem(eqs, t, vars, pars; name, systems, defaults = [input.u => area_int])
 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
+        p_b_int = p_b_int
+        d = d
+        x_int = x_int
+        Cd = Cd
+    end
+
+    systems = @named begin
+        port_a = HydraulicPort(; p_int = p_a_int)
+        port_b = HydraulicPort(; p_int = p_b_int)
+        flange = MechanicalPort()
+        valve = Valve(reversible; p_a_int, p_b_int,
+                      area_int = ParentScope(x_int) * 2π * ParentScope(d), Cd)
+    end
+
+    vars = @variables begin
+        x(t) = x_int
+        dx(t) = 0
+    end
+
+    eqs = [D(x) ~ dx
+           flange.v ~ dx
+           flange.f ~ 0 #TODO: model flow force
+           connect(valve.port_a, port_a)
+           connect(valve.port_b, port_b)
+           valve.input.u ~ x * 2π * d]
+
+    ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0])
+end
+
+@component function SpoolValve2Way(reversible = false; p_s_int, p_a_int, p_b_int, p_r_int,
+                                   m, g, x_int, Cd, d, name)
+    pars = @parameters begin
+        p_s_int = p_s_int
+        p_a_int = p_a_int
+        p_b_int = p_b_int
+        p_r_int = p_r_int
+
+        m = m
+        g = g
+
+        x_int = x_int
+
+        d = d
+
+        Cd = Cd
+    end
+
+    vars = []
+
+    systems = @named begin
+        vSA = SpoolValve(reversible; p_a_int = p_s_int, p_b_int = p_a_int, x_int, Cd, d)
+        vBR = SpoolValve(reversible; p_a_int = p_b_int, p_b_int = p_r_int, x_int, Cd, d)
+
+        port_s = HydraulicPort(; p_int = p_s_int)
+        port_a = HydraulicPort(; p_int = p_a_int)
+        port_b = HydraulicPort(; p_int = p_b_int)
+        port_r = HydraulicPort(; p_int = p_r_int)
+
+        mass = Mass(; m = m, g = g)
+
+        flange = MechanicalPort()
+    end
+
+    eqs = [connect(vSA.port_a, port_s)
+           connect(vSA.port_b, port_a)
+           connect(vBR.port_a, port_b)
+           connect(vBR.port_b, port_r)
+           connect(vSA.flange, vBR.flange, mass.flange, flange)]
+
+    ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0])
+end
+
+@component function Actuator(; p_a_int, p_b_int, area_a, area_b, length_a_int, length_b_int,
+                             m, g, x_int = 0, minimum_volume_a = 0, minimum_volume_b = 0,
+                             name)
+    pars = @parameters begin
+        p_a_int = p_a_int
+        p_b_int = p_b_int
+        area_a = area_a
+        area_b = area_b
+        x_int = x_int
+        length_a_int = length_a_int
+        length_b_int = length_b_int
+        minimum_volume_a = minimum_volume_a
+        minimum_volume_b = minimum_volume_b
+        m = m
+        g = g
+    end
+
+    vars = @variables begin
+        x(t) = x_int
+        dx(t) = 0
+    end
+
+    systems = @named begin
+        vol_a = DynamicVolume(; p_int = p_a_int, x_int = +x_int, area = area_a,
+                              dead_volume = length_a_int * area_a,
+                              minimum_volume = minimum_volume_a, direction = +1)
+        vol_b = DynamicVolume(; p_int = p_b_int, x_int = -x_int, area = area_b,
+                              dead_volume = length_b_int * area_b,
+                              minimum_volume = minimum_volume_b, direction = -1)
+        mass = Mass(; m, g)
+        port_a = HydraulicPort(; p_int = p_a_int)
+        port_b = HydraulicPort(; p_int = p_b_int)
+        flange = MechanicalPort()
+    end
+
+    eqs = [connect(vol_a.port, port_a)
+           connect(vol_b.port, port_b)
+           connect(vol_a.flange, vol_b.flange, mass.flange, flange)
+           x ~ vol_a.x
+           dx ~ vol_a.dx]
+
+    ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0])
+end
diff --git a/src/Hydraulic/IsothermalCompressible/utils.jl b/src/Hydraulic/IsothermalCompressible/utils.jl
index 4d112e27d..06b2b1f87 100644
--- a/src/Hydraulic/IsothermalCompressible/utils.jl
+++ b/src/Hydraulic/IsothermalCompressible/utils.jl
@@ -8,7 +8,7 @@ Connector port for hydraulic components.
 - `p_int`: [Pa] initial gauge pressure
 
 # States:
-- `p`: [Pa] gauge pressure
+- `p`: [Pa] gauge total pressure
 - `dm`: [kg/s] mass flow
 """
 @connector function HydraulicPort(; p_int, name)
@@ -16,6 +16,9 @@ Connector port for hydraulic components.
         ρ
         β
         μ
+        n
+        ρ_gas
+        p_gas
     end
 
     vars = @variables begin
@@ -29,20 +32,27 @@ end
 """
     HydraulicFluid(; density=997, bulk_modulus=2.09e9, viscosity=0.0010016, name)
 
-Fluid parameter setter for isothermal compressible fluid domain.  Defaults given for water at 20°C and 1atm reference pressure.
+Fluid parameter setter for isothermal compressible fluid domain.  Defaults given for water at 20°C and 0Pa gage (1atm absolute) reference pressure. Density is modeled using the Tait equation of state.  For pressures below the reference pressure, density is linearly interpolated to the gas state, this helps prevent pressures from going below the reference pressure.  
 
 # Parameters:
 
-- `ρ`: [kg/m^3] fluid density at reference pressure (set by `density` argument)
+- `ρ`: [kg/m^3] fluid density at 0Pa reference gage pressure (set by `density` argument)
 - `Β`: [Pa] fluid bulk modulus describing the compressibility (set by `bulk_modulus` argument)
 - `μ`: [Pa*s] or [kg/m-s] fluid dynamic viscosity  (set by `viscosity` argument)
+- `n`: density exponent
+- `ρ_gas`: [kg/m^3] density of fluid in gas state at reference gage pressure `p_gas` (set by `gas_density` argument)
+- `p_gas`: [Pa] reference pressure (set by `gas_pressure` argument)
 """
 @connector function HydraulicFluid(; density = 997, bulk_modulus = 2.09e9,
-                                   viscosity = 0.0010016, name)
+                                   viscosity = 0.0010016, gas_density = 0.0073955,
+                                   gas_pressure = -1000, n = 1, name)
     pars = @parameters begin
         ρ = density
         β = bulk_modulus
         μ = viscosity
+        n = n
+        ρ_gas = gas_density
+        p_gas = gas_pressure
     end
 
     vars = @variables begin dm(t), [connect = Flow] end
@@ -89,6 +99,8 @@ function friction_factor(dm, area, d_h, density, viscosity, shape_factor)
     return f
 end
 @register_symbolic friction_factor(dm, area, d_h, density, viscosity, shape_factor)
+Symbolics.derivative(::typeof(friction_factor), args, ::Val{1}) = 0
+Symbolics.derivative(::typeof(friction_factor), args, ::Val{4}) = 0
 
 function transition(x1, x2, y1, y2, x)
     if x < x1
@@ -102,7 +114,14 @@ function transition(x1, x2, y1, y2, x)
     end
 end
 
-density(port) = port.ρ
+density_ref(port) = port.ρ
+gas_density_ref(port) = port.ρ_gas
+gas_pressure_ref(port) = port.p_gas
 bulk_modulus(port) = port.β
 viscosity(port) = port.μ
-density(port, p) = density(port) * (1 + p / bulk_modulus(port))
+liquid_density(port) = density_ref(port) * (1 + port.p / bulk_modulus(port))
+function gas_density(port)
+    density_ref(port) -
+    port.p * (density_ref(port) - gas_density_ref(port)) / gas_pressure_ref(port)
+end
+full_density(port) = liquid_density(port) #ifelse( port.p > 0, liquid_density(port), gas_density(port) )
diff --git a/src/Mechanical/Translational/sources.jl b/src/Mechanical/Translational/sources.jl
index 73225e7c6..2e60580c6 100644
--- a/src/Mechanical/Translational/sources.jl
+++ b/src/Mechanical/Translational/sources.jl
@@ -2,12 +2,27 @@ using ModelingToolkitStandardLibrary.Blocks
 
 @component function Force(; name)
     @named flange = MechanicalPort()
-    @named f = RealInput()
+    @named input = RealInput()
 
     vars = pars = []
     eqs = [
-        flange.f ~ -f.u,
+        flange.f ~ -input.u,
     ]
     compose(ODESystem(eqs, t, vars, pars; name = name, defaults = [flange.v => 0]),
-            [flange, f])
+            [flange, input])
+end
+
+@component function Position(; x_int = 0, name)
+    @named flange = MechanicalPort()
+    @named input = RealInput()
+    pars = @parameters begin x_int = x_int end
+    vars = @variables begin
+        x(t) = x_int
+        dx(t) = 0
+    end
+    eqs = [x ~ input.u
+           D(x) ~ dx
+           flange.v ~ dx]
+    compose(ODESystem(eqs, t, vars, pars; name = name, defaults = [flange.v => 0]),
+            [flange, input])
 end
diff --git a/src/Mechanical/TranslationalPosition/sources.jl b/src/Mechanical/TranslationalPosition/sources.jl
index 02ad284df..0e11b4d20 100644
--- a/src/Mechanical/TranslationalPosition/sources.jl
+++ b/src/Mechanical/TranslationalPosition/sources.jl
@@ -6,7 +6,8 @@ Input signal acting as external force on a flange
 @component function Force(; name, use_support = false)
     @named partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
     @unpack flange = partial_element
-    @named f = RealInput() # Accelerating force acting at flange (= -flange.tau)
-    eqs = [flange.f ~ -f.u]
-    return extend(ODESystem(eqs, t, [], []; name = name, systems = [f]), partial_element)
+    @named input = RealInput() # Accelerating force acting at flange (= -flange.tau)
+    eqs = [flange.f ~ -input.u]
+    return extend(ODESystem(eqs, t, [], []; name = name, systems = [input]),
+                  partial_element)
 end
diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl
index 0e54a0f03..1c87f36d3 100644
--- a/test/Blocks/sources.jl
+++ b/test/Blocks/sources.jl
@@ -25,7 +25,7 @@ using OrdinaryDiffEq: ReturnCode.Success
 end
 
 @testset "TimeVaryingFunction" begin
-    f(t) = t^2
+    f(t) = t^2 + 1
 
     @named src = TimeVaryingFunction(f)
     @named int = Integrator()
@@ -41,7 +41,7 @@ end
     sol = solve(prob, Rodas4())
     @test sol.retcode == Success
     @test sol[src.output.u]≈f.(sol.t) atol=1e-3
-    @test sol[int.output.u][end]≈1 / 3 * 10^3 atol=1e-3 # closed-form solution to integral
+    @test sol[int.output.u][end]≈1 / 3 * 10^3 + 10 atol=1e-3 # closed-form solution to integral
 end
 
 @testset "Sine" begin
@@ -406,3 +406,30 @@ end
                                                              amplitude, damping, phase,
                                                              offset, start_time) atol=1e-3
 end
+
+@testset "Input" begin
+    dt = 4e-4
+    t_end = 10.0
+    time = 0:dt:t_end
+    x = @. time^2 + 1.0
+
+    @named src = Input(Float64)
+    @named int = Integrator()
+    @named iosys = ODESystem([
+                                 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
+
+    @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
+end
diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl
index f657b48f2..5d2075f05 100644
--- a/test/Hydraulic/isothermal_compressible.jl
+++ b/test/Hydraulic/isothermal_compressible.jl
@@ -3,149 +3,257 @@ import ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible as IC
 import ModelingToolkitStandardLibrary.Blocks as B
 import ModelingToolkitStandardLibrary.Mechanical.Translational as T
 
+using ModelingToolkitStandardLibrary.Blocks: Parameter
+
 @parameters t
 D = Differential(t)
 
-# ------------------------------------------------------------
-# Test fluid domain and Pipe
-function System(N; bulk_modulus, name)
-    pars = @parameters begin bulk_modulus = bulk_modulus end
-
-    systems = @named begin
-        fluid = IC.HydraulicFluid(; bulk_modulus)
-        stp = B.Step(; height = 10e5, offset = 0, start_time = 0.005, duration = Inf,
-                     smooth = true)
-        src = IC.InputSource(; p_int = 0)
-        vol = IC.FixedVolume(; p_int = 0, vol = 10.0)
-    end
+NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax = 4 // 10)
 
-    if N == 1
-        @named res = IC.PipeBase(; p_int = 0, area = 0.01, length = 500.0)
-    else
-        @named res = IC.Pipe(N; p_int = 0, area = 0.01, length = 500.0)
+@testset "Fluid Domain and Tube" begin
+    function System(N; bulk_modulus, name)
+        pars = @parameters begin bulk_modulus = bulk_modulus end
+
+        systems = @named begin
+            fluid = IC.HydraulicFluid(; bulk_modulus)
+            stp = B.Step(; height = 10e5, offset = 0, start_time = 0.005, duration = Inf,
+                         smooth = true)
+            src = IC.InputSource(; p_int = 0)
+            vol = IC.FixedVolume(; p_int = 0, vol = 10.0)
+        end
+
+        if N == 1
+            @named res = IC.TubeBase(; p_int = 0, area = 0.01, length = 500.0)
+        else
+            @named res = IC.Tube(N; p_int = 0, area = 0.01, length = 500.0)
+        end
+        push!(systems, res)
+
+        eqs = [connect(stp.output, src.input)
+               connect(fluid, src.port)
+               connect(src.port, res.port_a)
+               connect(res.port_b, vol.port)]
+
+        ODESystem(eqs, t, [], pars; name, systems)
     end
-    push!(systems, res)
 
-    eqs = [connect(stp.output, src.input)
-           connect(fluid, src.port)
-           connect(src.port, res.port_a)
-           connect(res.port_b, vol.port)]
+    @named sys1_2 = System(1; bulk_modulus = 2e9)
+    @named sys1_1 = System(1; bulk_modulus = 1e9)
+    @named sys5_1 = System(5; bulk_modulus = 1e9)
 
-    ODESystem(eqs, t, [], pars; name, systems)
-end
+    syss = structural_simplify.([sys1_2, sys1_1, sys5_1])
+    probs = [ODEProblem(sys, [], (0, 0.2)) for sys in syss] #ModelingToolkit.missing_variable_defaults(sys)
+    sols = [solve(prob, ImplicitEuler(nlsolve = NEWTON); initializealg = NoInit())
+            for prob in probs]
 
-@named sys1_2 = System(1; bulk_modulus = 2e9)
-@named sys1_1 = System(1; bulk_modulus = 1e9)
-@named sys5_1 = System(5; bulk_modulus = 1e9)
+    s1_2 = complete(sys1_2)
+    s1_1 = complete(sys1_1)
+    s5_1 = complete(sys5_1)
 
-NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax = 4 // 10)
+    # higher stiffness should compress more quickly and give a higher pressure
+    @test sols[1][s1_2.vol.port.p][end] > sols[2][s1_1.vol.port.p][end]
 
-syss = structural_simplify.([sys1_2, sys1_1, sys5_1])
-probs = [ODEProblem(sys, [], (0, 0.2)) for sys in syss];
-sols = [solve(prob, ImplicitEuler(nlsolve = NEWTON); initializealg = NoInit())
-        for prob in probs];
-
-s1_2 = complete(sys1_2)
-s1_1 = complete(sys1_1)
-s5_1 = complete(sys5_1)
-
-# higher stiffness should compress more quickly and give a higher pressure
-@test sols[1][s1_2.vol.port.p][end] > sols[2][s1_1.vol.port.p][end]
-
-# N=5 pipe is compressible, will pressurize more slowly
-@test sols[2][s1_1.vol.port.p][end] > sols[3][s5_1.vol.port.p][end]
-
-# ------------------------------------------------------------
-# Test the Valve 
-function System(; name)
-    pars = []
-
-    systems = @named begin
-        fluid = IC.HydraulicFluid()
-        sink = IC.Source(; p = 10e5)
-        vol = IC.FixedVolume(; vol = 0.1, p_int = 100e5)
-        valve = IC.Valve(; p_a_int = 10e5, p_b_int = 100e5, area_int = 0, Cd = 1e6)
-        ramp = B.Ramp(; height = 1, duration = 0.001, offset = 0, start_time = 0.001,
-                      smooth = true)
+    # N=5 pipe is compressible, will pressurize more slowly
+    @test sols[2][s1_1.vol.port.p][end] > sols[3][s5_1.vol.port.p][end]
+end
+
+@testset "Valve" begin
+    function System(; name)
+        pars = []
+
+        systems = @named begin
+            fluid = IC.HydraulicFluid()
+            sink = IC.Source(; p = 10e5)
+            vol = IC.FixedVolume(; vol = 0.1, p_int = 100e5)
+            valve = IC.Valve(; p_a_int = 10e5, p_b_int = 100e5, area_int = 0, Cd = 1e6)
+            ramp = B.Ramp(; height = 1, duration = 0.001, offset = 0, start_time = 0.001,
+                          smooth = true)
+        end
+
+        eqs = [connect(fluid, sink.port)
+               connect(sink.port, valve.port_a)
+               connect(valve.port_b, vol.port)
+               connect(valve.input, ramp.output)]
+
+        ODESystem(eqs, t, [], pars; name, systems)
     end
 
-    eqs = [connect(fluid, sink.port)
-           connect(sink.port, valve.port_a)
-           connect(valve.port_b, vol.port)
-           connect(valve.input, ramp.output)]
+    @named valve_system = System()
+    s = complete(valve_system)
+    sys = structural_simplify(valve_system)
+    prob = ODEProblem(sys, [], (0, 0.01))
+    sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4,
+                initializealg = NoInit())
 
-    ODESystem(eqs, t, [], pars; name, systems)
+    # the volume should discharge to 10bar
+    @test sol[s.vol.port.p][end] ≈ 10e5
 end
 
-@named valve_system = System()
-s = complete(valve_system)
-sys = structural_simplify(valve_system)
-prob = ODEProblem(sys, [], (0, 0.01))
-sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4,
-            initializealg = NoInit())
+@testset "DynamicVolume and minimum_volume feature" begin
+    function System(; name)
+        pars = []
 
-# the volume should discharge to 10bar
-@test sol[s.vol.port.p][end] ≈ 10e5
+        # DynamicVolume values
+        area = 0.01
+        length = 0.1
 
-# ------------------------------------------------------------
-# Test DynmaicVolume and minimum_volume feature
-function System(; name)
-    pars = []
+        systems = @named begin
+            fluid = IC.HydraulicFluid()
+            src1 = IC.InputSource(; p_int = 10e5)
+            src2 = IC.InputSource(; p_int = 10e5)
 
-    # DynamicVolume values
-    area = 0.01
-    length = 0.1
+            vol1 = IC.DynamicVolume(; p_int = 10e5, area, direction = +1,
+                                    dead_volume = 2e-4,
+                                    minimum_volume = 2e-4)
+            vol2 = IC.DynamicVolume(; p_int = 10e5, area, direction = -1,
+                                    dead_volume = 2e-4,
+                                    minimum_volume = 2e-4, x_int = length)
 
-    systems = @named begin
-        fluid = IC.HydraulicFluid()
-        src1 = IC.InputSource(; p_int = 10e5)
-        src2 = IC.InputSource(; p_int = 10e5)
+            mass = T.Mass(; m = 100)
 
-        vol1 = IC.DynamicVolume(; p_int = 10e5, area, direction = +1, dead_volume = 2e-4,
-                                minimum_volume = 2e-4)
-        vol2 = IC.DynamicVolume(; p_int = 10e5, area, direction = -1, dead_volume = 2e-4,
-                                minimum_volume = 2e-4, x_int = length)
+            sin1 = B.Sine(; frequency = 0.5, amplitude = +1e5, offset = 10e5)
+            sin2 = B.Sine(; frequency = 0.5, amplitude = -1e5, offset = 10e5)
+        end
 
-        mass = T.Mass(; m = 100)
+        eqs = [connect(fluid, src1.port, src2.port)
+               connect(src1.port, vol1.port)
+               connect(src2.port, vol2.port)
+               connect(vol1.flange, mass.flange, vol2.flange)
+               connect(src1.input, sin1.output)
+               connect(src2.input, sin2.output)]
 
-        sin1 = B.Sine(; frequency = 0.5, amplitude = +1e5, offset = 10e5)
-        sin2 = B.Sine(; frequency = 0.5, amplitude = -1e5, offset = 10e5)
+        ODESystem(eqs, t, [], pars; name, systems)
     end
 
-    eqs = [connect(fluid, src1.port, src2.port)
-           connect(src1.port, vol1.port)
-           connect(src2.port, vol2.port)
-           connect(vol1.flange, mass.flange, vol2.flange)
-           connect(src1.input, sin1.output)
-           connect(src2.input, sin2.output)]
-
-    ODESystem(eqs, t, [], pars; name, systems)
+    @named min_vol_system = System()
+    s = complete(min_vol_system)
+    sys = structural_simplify(min_vol_system)
+    prob = ODEProblem(sys, [], (0, 5.0))
+
+    sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4,
+                initializealg = NoInit())
+
+    # begin
+    #     fig = Figure()
+
+    #     ax = Axis(fig[1,1], ylabel="position [m]", xlabel="time [s]")
+    #     lines!(ax, sol.t, sol[s.vol1.x]; label="vol1")
+    #     lines!(ax, sol.t, sol[s.vol2.x]; label="vol2")
+
+    #     ax = Axis(fig[2,1], ylabel="pressure [bar]", xlabel="time [s]")
+    #     lines!(ax, sol.t, sol[s.vol1.port.p]/1e5; label="vol1")
+    #     lines!(ax, sol.t, sol[s.vol2.port.p]/1e5; label="vol2")
+    #     Legend(fig[2,2], ax)
+    #     fig    
+    # end
+
+    # volume/mass should stop moving at opposite ends
+    @test round(sol[s.vol1.x][1]; digits = 2) == 0.0
+    @test round(sol[s.vol2.x][1]; digits = 2) == 0.1
+    @test round(sol[s.vol1.x][end]; digits = 2) == 0.1
+    @test round(sol[s.vol2.x][end]; digits = 2) == 0.0
 end
 
-@named min_vol_system = System()
-s = complete(min_vol_system)
-sys = structural_simplify(min_vol_system)
-prob = ODEProblem(sys, [], (0, 5.0))
-
-sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt = 1e-4,
-            initializealg = NoInit())
-
-# begin
-#     fig = Figure()
-
-#     ax = Axis(fig[1,1], ylabel="position [m]", xlabel="time [s]")
-#     lines!(ax, sol.t, sol[s.vol1.x]; label="vol1")
-#     lines!(ax, sol.t, sol[s.vol2.x]; label="vol2")
-
-#     ax = Axis(fig[2,1], ylabel="pressure [bar]", xlabel="time [s]")
-#     lines!(ax, sol.t, sol[s.vol1.port.p]/1e5; label="vol1")
-#     lines!(ax, sol.t, sol[s.vol2.port.p]/1e5; label="vol2")
-#     Legend(fig[2,2], ax)
-#     fig    
-# end
-
-# volume/mass should stop moving at opposite ends
-@test round(sol[s.vol1.x][1]; digits = 2) == 0.0
-@test round(sol[s.vol2.x][1]; digits = 2) == 0.1
-@test round(sol[s.vol1.x][end]; digits = 2) == 0.1
-@test round(sol[s.vol2.x][end]; digits = 2) == 0.0
+@testset "Actuator System" begin
+    function System(use_input, f; name)
+        @parameters t
+
+        pars = @parameters begin
+            p_s = 285e5
+            p_r = 5e5
+
+            A_1 = 908e-4
+            A_2 = 360e-4
+
+            p_1 = 17e5
+            p_2 = p_1 * A_1 / A_2
+
+            l_1 = 0.7
+            l_2 = 2.7
+            m_f = 276
+            g = 0
+            x_f_int = 0
+
+            d = 230e-3
+
+            Cd = 70e5 * 2 * 876 / (876 * 50 / (25e-3 * 2π * 230e-3))^2   #50_000 lpm
+
+            m_piston = 880
+            m_body = 1500
+        end
+
+        vars = @variables begin ddx(t) = 0 end
+
+        systems = @named begin
+            src = IC.Source(; p = p_s)
+            valve = IC.SpoolValve2Way(; p_s_int = p_s, p_a_int = p_1, p_b_int = p_2,
+                                      p_r_int = p_r, g, m = m_f, x_int = x_f_int, d, Cd)
+            piston = IC.Actuator(; p_a_int = p_1, p_b_int = p_2, area_a = A_1, area_b = A_2,
+                                 length_a_int = l_1, length_b_int = l_2, m = m_piston,
+                                 g = 0,
+                                 x_int = 0, minimum_volume_a = 0, minimum_volume_b = 0)
+            body = T.Mass(; m = m_body)
+            pipe = IC.Tube(5; p_int = p_2, area = A_2, length = 2.0)
+            snk = IC.Source(; p = p_r)
+            pos = T.Position(; x_int = x_f_int)
+
+            m1 = IC.FlowDivider(; p_int = p_2, n = 3)
+            m2 = IC.FlowDivider(; p_int = p_2, n = 3)
+
+            fluid = IC.HydraulicFluid()
+        end
+
+        if use_input
+            @named input = B.Input(Float64)
+        else
+            @named input = B.TimeVaryingFunction(f; t)
+        end
+
+        push!(systems, input)
+
+        eqs = [connect(input.output, pos.input)
+               connect(valve.flange, pos.flange)
+               connect(valve.port_a, piston.port_a)
+               connect(piston.flange, body.flange)
+               connect(piston.port_b, m1.port_a)
+               connect(m1.port_b, pipe.port_b)
+               connect(pipe.port_a, m2.port_b)
+               connect(m2.port_a, valve.port_b)
+               connect(src.port, valve.port_s)
+               connect(snk.port, valve.port_r)
+               connect(fluid, src.port, snk.port)
+               D(body.v) ~ ddx]
+
+        ODESystem(eqs, t, vars, pars; name, systems)
+    end
+
+    @named system = System(true, nothing)
+
+    sys = structural_simplify(system)
+    defs = ModelingToolkit.defaults(sys)
+    s = complete(system)
+    prob = ODEProblem(sys, [], (0, 1.0); tofloat = false, jac = true)
+
+    # check the fluid domain
+    @test Symbol(defs[s.src.port.ρ]) == Symbol(s.fluid.ρ)
+    @test Symbol(defs[s.valve.port_s.ρ]) == Symbol(s.fluid.ρ)
+    @test Symbol(defs[s.valve.port_a.ρ]) == Symbol(s.fluid.ρ)
+    @test Symbol(defs[s.valve.port_b.ρ]) == Symbol(s.fluid.ρ)
+    @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.055
+    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)
+
+    p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false))
+    prob = remake(prob; p, tspan = (0, time[end]))
+    sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt,
+                initializealg = NoInit())
+
+    @test sol[sys.ddx][1] == 0.0
+    @test maximum(sol[sys.ddx]) > 500
+    @test sol[sys.ddx][end] < 100
+end
diff --git a/test/Mechanical/translational.jl b/test/Mechanical/translational.jl
index b9729752a..b2e3bb3e8 100644
--- a/test/Mechanical/translational.jl
+++ b/test/Mechanical/translational.jl
@@ -65,7 +65,7 @@ end
     @named source = Sine(frequency = 3, amplitude = 2)
 
     function simplify_and_solve(damping, spring, body, ground, f, source)
-        eqs = [connect(f.f, source.output)
+        eqs = [connect(f.input, source.output)
                connect(f.flange, body.flange)
                connect(spring.flange_a, body.flange, damping.flange_a)
                connect(spring.flange_b, damping.flange_b, ground.flange)]