fix test

Author: ValentinKaisermayer

src/Hydraulic/components.jl

@@ -12,10 +12,8 @@ Infinite reservoir; set BC for hydraulic system
 function Reservoir(; name, p)
     @parameters p = p
     @named a = HydraulicPort()
-    eqs = [
-        a.p ~ p
-    ]
-    return ODESystem(eqs, t, [], [p]; systems=[a], name=name)
+    eqs = [a.p ~ p]
+    return ODESystem(eqs, t, [], [p]; systems = [a], name = name)
 end
 
 """
@@ -32,13 +30,13 @@ end
 - `Cd`: [] Discharge coefficient
 - `Re_crit=150`: [] Critical Reynolds number
 """
-function LocalRestriction(; name, A, Cd, Re_crit=150)
+function LocalRestriction(; name, A, Cd, Re_crit = 150)
     @named a = HydraulicPort()
     @named b = HydraulicPort()
 
     pars = @parameters begin
-        A=A 
-        Cd=Cd 
+        A = A
+        Cd = Cd
         p_atm# , [scope=:parent]
         nu_atm# , [scope=:parent]
         beta_atm# , [scope=:parent]
@@ -55,15 +53,15 @@ function LocalRestriction(; name, A, Cd, Re_crit=150)
     rho = calc_density((b.p + a.p) / 2, rho_atm, p_atm, beta_atm)
 
     Δp = b.p - a.p
-    p_cr = pi/4 * rho / (2 * A) * (Re_crit * nu_atm / Cd)^2
+    p_cr = pi / 4 * rho / (2 * A) * (Re_crit * nu_atm / Cd)^2
 
     eqs = [
         # Momentum balance
         b.m_flow ~ Cd * A * sqrt(2 * rho) * regRoot(Δp, p_cr)
         # Mass balance
         a.m_flow + b.m_flow ~ 0
     ]
-    return ODESystem(eqs, t, [], pars; name=name, systems=[a, b])
+    return ODESystem(eqs, t, [], pars; name = name, systems = [a, b])
 end
 
 """
@@ -81,11 +79,11 @@ Constant volume chamber.
 - `V`: [`m^3`] Constant volume
 - `p_start=0.0`: [`Pa`] Initial pressure inside the volume
 """
-function ConstantVolume(; name, V, p_start=0.0)
+function ConstantVolume(; name, V, p_start = 0.0)
     @named a = HydraulicPort()
-    @variables p(t)=p_start # pressure inside the volume
+    @variables p(t) = p_start # pressure inside the volume
     pars = @parameters begin
-        V=V
+        V = V
         p_atm# , [scope=:parent]
         nu_atm# , [scope=:parent]
         beta_atm# , [scope=:parent]
@@ -100,12 +98,12 @@ function ConstantVolume(; name, V, p_start=0.0)
     pars = [V, p_atm, nu_atm, beta_atm, rho_atm]
 
     rho = calc_density(p, rho_atm, p_atm, beta_atm)
-    
+
     eqs = [
         # Mass conservation
         D(p) ~ 1 / (Symbolics.derivative(rho, p) * V) * a.m_flow
         a.p ~ p # no pressure loss from inlet to volume
     ]
 
-    return ODESystem(eqs, t, [p], pars; name=name, systems=[a])
-end
+    return ODESystem(eqs, t, [p], pars; name = name, systems = [a])
+end

src/Hydraulic/sources.jl

@@ -18,7 +18,7 @@ function PressureSource(; name, dp)
         a.p ~ b.p + dp
         a.m_flow + b.m_flow ~ 0
     ]
-    return ODESystem(eqs, t, [], [dp]; systems=[a, b], name=name)
+    return ODESystem(eqs, t, [], [dp]; systems = [a, b], name = name)
 end
 
 """
@@ -41,5 +41,5 @@ function MassFlowRateSource(; name, m_flow)
         a.m_flow + b.m_flow ~ 0
         a.m_flow ~ m_flow
     ]
-    return ODESystem(eqs, t, [], [m_flow]; systems=[a, b], name=name)
-end
+    return ODESystem(eqs, t, [], [m_flow]; systems = [a, b], name = name)
+end

src/Hydraulic/utils.jl

@@ -1,7 +1,7 @@
-@connector function HydraulicPort(; name, p_start=1e5, m_flow_start=0.0)
+@connector function HydraulicPort(; name, p_start = 1e5, m_flow_start = 0.0)
     @variables p(t) = p_start
     @variables m_flow(t) = m_flow_start [connect = Flow]
-    return ODESystem(Equation[], t, [p, m_flow], []; name=name)
+    return ODESystem(Equation[], t, [p, m_flow], []; name = name)
 end
 @doc """
     HydraulicPort(; name, p_start=1e5, m_flow_start=0.0)
@@ -18,10 +18,10 @@ Port for hydraulic systems.
 """ HydraulicPort
 
 # smooth approximation of x * abs(x)
-regAbs(x, delta=0.01) = x * sqrt(x^2 + delta^2)
+regAbs(x, delta = 0.01) = x * sqrt(x^2 + delta^2)
 
 # Smoothed version of sign(x)*sqrt(abs(x)), which has a finite derivative at x=0
-regRoot(x, delta=0.01) = x / (x^2 + delta^2)^0.25
+regRoot(x, delta = 0.01) = x / (x^2 + delta^2)^0.25
 
 """
     FluidProperties(;name, p_atm=101325.0, nu_atm=1.0034e-6, beta_atm=2.1791e9, rho_atm=998.21)
@@ -40,9 +40,15 @@ model = extend(model, fluid_props)
 - `beta_atm = 2.1791e9`: [Pa] Isothermal bulk modulus at atmospheric pressure
 - `rho_atm = 998.21`: [kg/m^3] Liquid density at atmospheric pressure
 """
-function FluidProperties(;name, p_atm=101325.0, nu_atm=1.0034e-6, beta_atm=2.1791e9, rho_atm=998.21)
-    pars = @parameters p_atm=p_atm nu_atm=nu_atm beta_atm=beta_atm rho_atm=rho_atm
-    ODESystem(Equation[], t, [], pars; name=name)
+function FluidProperties(;
+    name,
+    p_atm = 101325.0,
+    nu_atm = 1.0034e-6,
+    beta_atm = 2.1791e9,
+    rho_atm = 998.21,
+)
+    pars = @parameters p_atm = p_atm nu_atm = nu_atm beta_atm = beta_atm rho_atm = rho_atm
+    ODESystem(Equation[], t, [], pars; name = name)
 end
 
 """
@@ -65,4 +71,4 @@ Assumptions:
 """
 function calc_density(p, rho_atm, p_atm, beta_atm)
     return rho_atm * exp((p - p_atm) / beta_atm)
-end 
+end

test/Hydraulic/hydraulic.jl

@@ -6,29 +6,33 @@ using ModelingToolkit, OrdinaryDiffEq, Test
 D = Differential(t)
 
 @testset "restriction and volume (shows dynamic compressibility)" begin
-    @named reservoir_left = Reservoir(; p=2e5)
-    @named valve = LocalRestriction(; A=1e-4, Cd=0.64)
-    @named volume = ConstantVolume(; V=0.001)
+    @named reservoir_left = Reservoir(; p = 2e5)
+    @named valve = LocalRestriction(; A = 1e-4, Cd = 0.64)
+    @named volume = ConstantVolume(; V = 0.001)
     connections = [
         connect(reservoir_left.a, valve.a)
         connect(valve.b, volume.a)
     ]
-    @named model = 
-        ODESystem(connections, t; systems=[reservoir_left, valve, volume])
-    
+    @named model = ODESystem(connections, t; systems = [reservoir_left, valve, volume])
+
     # Water
-    @named fluid_props = 
-        FluidProperties(;p_atm=101325.0, nu_atm=1.0034e-6, beta_atm=2.1791e9, rho_atm=998.21)
+    @named fluid_props = FluidProperties(;
+        p_atm = 101325.0,
+        nu_atm = 1.0034e-6,
+        beta_atm = 2.1791e9,
+        rho_atm = 998.21,
+    )
     model = extend(model, fluid_props)
 
     sys = structural_simplify(model)
     prob = ODEProblem(sys, [volume.p => 1e5], (0.0, 0.5e-3))
-    sol = solve(prob, Rodas4())
+    sol = solve(prob, Tsit5(), reltol = 1e-6)
 
     @test sol.retcode == :Success
-    @test sol[volume.p][end] ≈ 2e5
+    @test sol[volume.p][end] ≈ 2e5 atol = 1
+    @test sol[volume.a.m_flow][end] ≈ 0 atol = 1e-2
 
     # p1=plot(sol; vars=[volume.p / 1e6], ylabel="p in MPa", label="Pressure of Volume")
-    # p2=plot(sol; vars=[volume.a.m_flow], ylabel="m_flow in kg/s", label="Massflow rate into the volume")
+    # p2=plot(sol; vars=[volume.a.m_flow], ylabel="m_flow in kg/s", label="Mass Flow Rate into the volume")
     # plot(p1, p2, layout=(2,1))
 end