diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl
index 56228c8a..ce9b39ab 100644
--- a/src/Hydraulic/IsothermalCompressible/components.jl
+++ b/src/Hydraulic/IsothermalCompressible/components.jl
@@ -192,20 +192,17 @@ Constant length internal flow model discretized by `N` (`FixedVolume`: `N`, `Tub
         push!(pipe_bases, x)
     end
 
-    
-    eqs = [
-        connect(pipe_bases[1].port_a, port_a)
-        connect(pipe_bases[end].port_b, port_b)
-    ]
+    eqs = [connect(pipe_bases[1].port_a, port_a)
+           connect(pipe_bases[end].port_b, port_b)]
 
     volumes = []
-    for i in 1:(N-1)
+    for i in 1:(N - 1)
         x = FixedVolume(; name = Symbol("v$i"),
-            vol = ParentScope(area) * ParentScope(length) / (N-1),
-            p_int = ParentScope(p_int) )
+            vol = ParentScope(area) * ParentScope(length) / (N - 1),
+            p_int = ParentScope(p_int))
         push!(volumes, x)
         push!(eqs,
-            connect(x.port, pipe_bases[i].port_b, pipe_bases[i+1].port_a))
+            connect(x.port, pipe_bases[i].port_b, pipe_bases[i + 1].port_a))
     end
 
     return ODESystem(eqs, t, vars, pars; name, systems = [ports; pipe_bases; volumes])
@@ -268,8 +265,8 @@ end
     end
 
     vars = @variables begin
-        area(t), [guess=0]
-        y(t), [guess=0]
+        area(t), [guess = 0]
+        y(t), [guess = 0]
     end
 
     # let
@@ -361,7 +358,7 @@ end
         dx(t), [guess = 0]
         rho(t), [guess = liquid_density(port)]
         m(t), [guess = 0]
-        vol(t) 
+        vol(t)
     end
 
     # let
@@ -403,16 +400,15 @@ Fixed fluid volume.
 
     vars = @variables begin
         rho(t), [guess = liquid_density(port)]
-        m(t), [guess=vol*liquid_density(port)]
-        p(t)=p_int
+        m(t), [guess = vol * liquid_density(port)]
+        p(t) = p_int
     end
 
     # let
     dm = port.dm
-    
 
     eqs = [D(m) ~ dm
-        #    rho ~ full_density(port, p)
+           #    rho ~ full_density(port, p)
            p ~ full_pressure(port, rho) # see https://github.com/SciML/OrdinaryDiffEq.jl/issues/2561
            p ~ port.p
            m ~ rho * vol]
@@ -464,21 +460,20 @@ See also [`FixedVolume`](@ref), [`DynamicVolume`](@ref)
         area,
         direction = +1,
         x_int,
-         name)
-
+        name)
     pars = @parameters begin
         area = area
         x_int = x_int
     end
 
     vars = @variables begin
-        x(t)=x_int
-        dx(t), [guess=0]
-        p(t), [guess=0]
-        f(t), [guess=0]
-        rho(t), [guess=0]
-        m(t), [guess=0]
-        dm(t), [guess=0]
+        x(t) = x_int
+        dx(t), [guess = 0]
+        p(t), [guess = 0]
+        f(t), [guess = 0]
+        rho(t), [guess = 0]
+        m(t), [guess = 0]
+        dm(t), [guess = 0]
     end
 
     systems = @named begin
@@ -504,7 +499,6 @@ See also [`FixedVolume`](@ref), [`DynamicVolume`](@ref)
            # physics
            # rho ~ full_density(port, p)
            p ~ full_pressure(port, rho) # see https://github.com/SciML/OrdinaryDiffEq.jl/issues/2561
-
            f ~ p * area
            m ~ rho * x * area]
 
@@ -574,10 +568,7 @@ dm ────►               │  │ area
         minimum_area = 0,
 
         # Damping
-        d = 0,
-
-        name)
-
+        d = 0, name)
     @assert (direction == +1)||(direction == -1) "direction argument must be +/-1, found $direction"
 
     #TODO: How to set an assert effective_length >= length ??
@@ -602,9 +593,9 @@ dm ────►               │  │ area
     end
 
     vars = @variables begin
-        x(t)=x_int 
-        vol(t), [guess=x_int * area]
-    end 
+        x(t) = x_int
+        vol(t), [guess = x_int * area]
+    end
 
     systems = @named begin
         port = HydraulicPort(;)
@@ -617,7 +608,7 @@ dm ────►               │  │ area
             area,
             dead_volume = area * x_int,
             p_int,
-            x_int=0) 
+            x_int = 0)
     end
 
     ratio = (x - x_min) / (x_damp - x_min)
@@ -637,9 +628,7 @@ dm ────►               │  │ area
            connect(port, damper.port_b)
            connect(moving_volume.port, damper.port_a)
            dx ~ flange.v * direction
-           p * area - dx*d ~ -flange.f * direction
-           ]
-
+           p * area - dx * d ~ -flange.f * direction]
 
     return ODESystem(eqs, t, vars, pars; name, systems)
 end
@@ -678,8 +667,8 @@ See [`Valve`](@ref) for more information.
     end
 
     vars = @variables begin
-        x(t)=x_int
-        dx(t), [guess=0]
+        x(t) = x_int
+        dx(t), [guess = 0]
     end
 
     eqs = [D(x) ~ dx
@@ -751,10 +740,10 @@ See [`SpoolValve`](@ref) for more information.
            connect(vBR.port_a, port_b)
            connect(vBR.port_b, port_r)
            connect(vSA.flange, vBR.flange, mass.flange, flange)]
-    
+
     initialization_eqs = [
         mass.s ~ x_int
-        # mass.v ~ dx_int
+    # mass.v ~ dx_int
     ]
 
     ODESystem(eqs, t, vars, pars; name, systems, initialization_eqs)
@@ -848,7 +837,7 @@ Actuator made of two DynamicVolumes connected in opposite direction with body ma
         damping_volume_b = minimum_volume_b,
         Cd = 1e4,
         Cd_reverse = Cd,
-        d=0,
+        d = 0,
         p_a_int,
         p_b_int,
         name)
@@ -963,14 +952,14 @@ dm ────►  effective area
 
 @mtkmodel Orifice begin
     @parameters begin
-        orifice_area=0.00094
-        Cd=0.6 # TODO Cd here is defined differently from Valve(). 
+        orifice_area = 0.00094
+        Cd = 0.6 # TODO Cd here is defined differently from Valve(). 
         # Here it follows the form Effective Orifice Area = Cd x Physical Orifice Area
         # The Valve component should be updated too.
     end
     @components begin
-        area = Constant(k=orifice_area)
-        valve = Valve(Cd= 1 / (Cd * Cd))
+        area = Constant(k = orifice_area)
+        valve = Valve(Cd = 1 / (Cd * Cd))
         port_a = HydraulicPort()
         port_b = HydraulicPort()
     end
diff --git a/src/Hydraulic/IsothermalCompressible/utils.jl b/src/Hydraulic/IsothermalCompressible/utils.jl
index 2ba04755..1bad646e 100644
--- a/src/Hydraulic/IsothermalCompressible/utils.jl
+++ b/src/Hydraulic/IsothermalCompressible/utils.jl
@@ -1,5 +1,7 @@
 # regPow(x, a, delta = 0.01) = x * (x * x + delta * delta)^((a - 1) / 2);
-regPow(x, a, delta = 0.01) = ifelse(abs(x/delta) >= 1, sign(x)*abs(x/delta)^a * delta^a, (delta^a*x)/delta)
+function regPow(x, a, delta = 0.01)
+    ifelse(abs(x / delta) >= 1, sign(x) * abs(x / delta)^a * delta^a, (delta^a * x) / delta)
+end
 regRoot(x, delta = 0.01) = regPow(x, 0.5, delta)
 
 """
@@ -73,7 +75,6 @@ Fluid parameter setter for isothermal compressible fluid domain.  Defaults given
     ODESystem(eqs, t, vars, pars; name)
 end
 
-
 function transition(x1, x2, y1, y2, x)
     u = (x - x1) / (x2 - x1)
     blend = u^2 * (3 - 2 * u)
@@ -142,12 +143,9 @@ liquid_density(port) = liquid_density(port, port.p)
 # (p/beta + 1)*rho_0 = rho
 
 function liquid_pressure(port, rho)
-
-    (rho/density_ref(port) - 1)*bulk_modulus(port)
-
+    (rho / density_ref(port) - 1) * bulk_modulus(port)
 end
 
-
 function gas_density(port, p)
     slope = (density_ref(port) - gas_density_ref(port)) / (0 - gas_pressure_ref(port))
     b = density_ref(port)
@@ -156,8 +154,7 @@ function gas_density(port, p)
 end
 
 function gas_pressure(port, rho)
-
-    slope =  (0 - gas_pressure_ref(port)) / (density_ref(port) - gas_density_ref(port))
+    slope = (0 - gas_pressure_ref(port)) / (density_ref(port) - gas_density_ref(port))
     b = 0
 
     return b + rho * slope
@@ -172,7 +169,8 @@ full_density(port) = full_density(port, port.p)
 
 function full_pressure(port, rho)
     ifelse(port.let_gas == 1,
-        ifelse( rho >= density_ref(port), liquid_pressure(port, rho), gas_pressure(port, rho)),
+        ifelse(
+            rho >= density_ref(port), liquid_pressure(port, rho), gas_pressure(port, rho)),
         liquid_pressure(port, rho)
     )
-end
\ No newline at end of file
+end
diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl
index e4bda065..06c92725 100644
--- a/test/Hydraulic/isothermal_compressible.jl
+++ b/test/Hydraulic/isothermal_compressible.jl
@@ -18,20 +18,17 @@ NEWTON = NLNewton(
 
         systems = @named begin
             fluid = IC.HydraulicFluid(; bulk_modulus)
-            stp = B.Step(; height = 10e5, offset = 0, start_time = 0.005, duration = Inf, smooth = true)
+            stp = B.Step(; height = 10e5, offset = 0, start_time = 0.005,
+                duration = Inf, smooth = true)
             src = IC.Pressure()
             vol = IC.FixedVolume(; vol = 10.0, p_int)
             res = IC.Tube(N; area = 0.01, length = 50.0, p_int)
         end
 
-        eqs = [
-               connect(stp.output, src.p)
+        eqs = [connect(stp.output, src.p)
                connect(fluid, src.port)
-            
                connect(src.port, res.port_a)
-               connect(res.port_b, vol.port)
-               
-               ]
+               connect(res.port_b, vol.port)]
 
         ODESystem(eqs, t, [], pars; name, systems)
     end
@@ -40,8 +37,6 @@ NEWTON = NLNewton(
     @mtkbuild s1_2 = System(1; bulk_modulus = 2e9)
     @mtkbuild s5_1 = System(5; bulk_modulus = 1e9)
 
-    
-    
     p1_1 = ODEProblem(s1_1, [], (0, 0.05))
     p1_2 = ODEProblem(s1_2, [], (0, 0.05))
     p5_1 = ODEProblem(s5_1, [], (0, 0.05))
@@ -57,27 +52,25 @@ NEWTON = NLNewton(
     # lines!(ax, tm, sol1_2.(tm; idxs=s1_1.vol.port.p)); fig
     # lines!(ax, tm, sol5_1.(tm; idxs=s5_1.vol.port.p)); fig
     # fig
-    
+
     # higher stiffness should compress more quickly and give a higher pressure
     @test sol1_2[s1_2.vol.port.p][end] > sol1_1[s1_1.vol.port.p][end]
 
     # N=5 pipe is compressible, will pressurize more slowly
     @test sol1_1[s1_1.vol.port.p][end] > sol5_1[s5_1.vol.port.p][end]
-
-    
 end
 
 @testset "Valve" begin
-
     function System(; name)
         pars = []
 
         systems = @named begin
             fluid = IC.HydraulicFluid()
             sink = IC.FixedPressure(; p = 10e5)
-            vol = IC.FixedVolume(; vol = 5, p_int=1e5)
+            vol = IC.FixedVolume(; vol = 5, p_int = 1e5)
             valve = IC.Valve(; Cd = 1e5, minimum_area = 0)
-            ramp = B.Ramp(; height = 0.1, duration = 0.1, offset = 0, start_time = 0.1, smooth = true)
+            ramp = B.Ramp(;
+                height = 0.1, duration = 0.1, offset = 0, start_time = 0.1, smooth = true)
         end
 
         eqs = [connect(fluid, sink.port)
@@ -89,9 +82,9 @@ end
     end
 
     @named valve_system = System()
-    sys = structural_simplify(valve_system)  
+    sys = structural_simplify(valve_system)
     prob = ODEProblem(sys, [], (0, 1))
-    sol = solve(prob, Rodas5P(); abstol=1e-6, reltol=1e-9)
+    sol = solve(prob, Rodas5P(); abstol = 1e-6, reltol = 1e-9)
     s = complete(valve_system)
 
     # the volume should discharge to 10bar
@@ -105,7 +98,7 @@ end
 end
 
 @testset "DynamicVolume and minimum_volume feature" begin # Need help here
-    function System(; name, area = 0.01, length = 0.1, damping_volume=length*area*0.1)
+    function System(; name, area = 0.01, length = 0.1, damping_volume = length * area * 0.1)
         pars = []
 
         # DynamicVolume values
@@ -149,18 +142,16 @@ end
                connect(src1.p, sin1.output)
                connect(src2.p, sin2.output)]
 
-        initialization_eqs = [
-            mass.s ~ 0.0
-            mass.v ~ 0.0
-        ]
+        initialization_eqs = [mass.s ~ 0.0
+                              mass.v ~ 0.0]
 
         ODESystem(eqs, t, [], pars; name, systems, initialization_eqs)
     end
 
     @mtkbuild sys = System()
     prob = ODEProblem(sys, [], (0, 5))
-    sol = solve(prob, Rodas5P(); abstol=1e-6, reltol=1e-9)
-        
+    sol = solve(prob, Rodas5P(); abstol = 1e-6, reltol = 1e-9)
+
     # begin
     #     fig = Figure()
 
@@ -182,21 +173,20 @@ end
     # end
 
     # volume/mass should stop moving at opposite ends
-    @test sol(0; idxs=sys.vol1.x) == 0.1
-    @test sol(0; idxs=sys.vol2.x) == 0.1
+    @test sol(0; idxs = sys.vol1.x) == 0.1
+    @test sol(0; idxs = sys.vol2.x) == 0.1
 
-    @test round(sol(1; idxs=sys.vol1.x); digits=2) == 0.19
-    @test round(sol(1; idxs=sys.vol2.x); digits=2) == 0.01
+    @test round(sol(1; idxs = sys.vol1.x); digits = 2) == 0.19
+    @test round(sol(1; idxs = sys.vol2.x); digits = 2) == 0.01
 
-    @test round(sol(2; idxs=sys.vol1.x); digits=2) == 0.01
-    @test round(sol(2; idxs=sys.vol2.x); digits=2) == 0.19
+    @test round(sol(2; idxs = sys.vol1.x); digits = 2) == 0.01
+    @test round(sol(2; idxs = sys.vol2.x); digits = 2) == 0.19
 
-    @test round(sol(3; idxs=sys.vol1.x); digits=2) == 0.19
-    @test round(sol(3; idxs=sys.vol2.x); digits=2) == 0.01
-    
-    @test round(sol(4; idxs=sys.vol1.x); digits=2) == 0.01
-    @test round(sol(4; idxs=sys.vol2.x); digits=2) == 0.19
+    @test round(sol(3; idxs = sys.vol1.x); digits = 2) == 0.19
+    @test round(sol(3; idxs = sys.vol2.x); digits = 2) == 0.01
 
+    @test round(sol(4; idxs = sys.vol1.x); digits = 2) == 0.01
+    @test round(sol(4; idxs = sys.vol2.x); digits = 2) == 0.19
 end
 
 @testset "Actuator System" begin
@@ -258,7 +248,7 @@ end
             @named input = B.SampledData(Float64)
         else
             #@named input = B.TimeVaryingFunction(f)
-            @named input = B.Constant(k=0)
+            @named input = B.Constant(k = 0)
         end
 
         push!(systems, input)
@@ -266,39 +256,37 @@ end
         eqs = [connect(input.output, pos.s)
                connect(valve.flange, pos.flange)
                connect(valve.port_a, piston.port_a)
-            #    connect(piston.flange, body.flange)
+               #    connect(piston.flange, body.flange)
 
                connect(piston.port_b, valve.port_b)
-               
-            #    connect(piston.port_b, pipe.port_b)
-            # #    connect(piston.port_b, m1.port_a)
-            # #    connect(m1.port_b, pipe.port_b)
-            
-            # connect(pipe.port_a, valve.port_b)
-            # #    connect(pipe.port_a, m2.port_b)
-            # #    connect(m2.port_a, valve.port_b)
+
+               #    connect(piston.port_b, pipe.port_b)
+               # #    connect(piston.port_b, m1.port_a)
+               # #    connect(m1.port_b, pipe.port_b)
+
+               # connect(pipe.port_a, valve.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(piston.mass.v) ~ ddx
-               ]
-        
+               D(piston.mass.v) ~ ddx]
+
         initialization_eqs = [
-                # body.s ~ 0
-            ]
+        # body.s ~ 0
+        ]
 
         ODESystem(eqs, t, vars, pars; name, systems, initialization_eqs)
     end
 
     @mtkbuild initsys = System(false)
-    
-    initprob = ODEProblem(initsys, [], (0,0))
+
+    initprob = ODEProblem(initsys, [], (0, 0))
     initsol = solve(initprob, Rodas5P())
 
     @mtkbuild sys = System(true)
 
-
     dt = 1e-4
     time = 0:dt:0.1
 
@@ -310,7 +298,7 @@ end
 
     # NOTE: bypassing initialization system: https://github.com/SciML/ModelingToolkit.jl/issues/3312
     prob = ODEProblem(sys, initsol[1], (0, 0.1); build_initializeprob = false)
-       
+
     #TODO: Implement proper initialization system after issue is resolved
     #TODO: How to bring the body back and not have an overdetermined system?
 
@@ -322,16 +310,13 @@ end
     @test Symbol(defs[sys.valve.port_r.ρ]) == Symbol(sys.fluid.ρ)
     @test Symbol(defs[sys.snk.port.ρ]) == Symbol(sys.fluid.ρ)
 
-    @time sol = solve(prob, Rodas5P(); initializealg = NoInit());
-        
+    @time sol = solve(prob, Rodas5P(); initializealg = NoInit())
+
     @test sol[sys.ddx][1] == 0.0
     @test maximum(sol[sys.ddx]) > 200
     @test sol[sys.piston.x][end] > 0.6
-
-
 end
 
-
 @testset "Prevent Negative Pressure" begin
     @component function System(; name)
         pars = @parameters let_gas = 1
@@ -339,7 +324,7 @@ end
         systems = @named begin
             fluid = IC.HydraulicFluid(; let_gas)
             vol = IC.DynamicVolume(; area = 0.001, x_int = 0.05,
-                x_max = 0.1, x_damp = 0.02, x_min = 0.01, direction = +1, p_int=100e5)
+                x_max = 0.1, x_damp = 0.02, x_min = 0.01, direction = +1, p_int = 100e5)
             mass = T.Mass(; m = 100, g = -9.807) # s = 0.05
             cap = IC.Cap()
         end
@@ -347,23 +332,20 @@ end
         eqs = [connect(fluid, cap.port, vol.port)
                connect(vol.flange, mass.flange)]
 
-        initialization_eqs = [
-            mass.s ~ 0.05
-            mass.v ~ 0
-        ]
-            
+        initialization_eqs = [mass.s ~ 0.05
+                              mass.v ~ 0]
 
         return ODESystem(eqs, t, [], pars; name, systems, initialization_eqs)
     end
 
     @mtkbuild sys = System()
-    
+
     prob1 = ODEProblem(sys, [], (0, 0.05))
     # prob1 = remake(prob1; u0 = BigFloat.(prob1.u0))
     prob2 = ODEProblem(sys, [], (0, 0.05), [sys.let_gas => 0])
 
     # @time sol1 = solve(prob1, Rodas5P(); abstol=1e-9, reltol=1e-9) #BUG: Using BigFloat gives... ERROR: MethodError: no method matching getindex(::Missing, ::Int64)
-    @time sol1 = solve(prob1, Rodas5P(); adaptive=false, dt=1e-6) #TODO: fix BigFloat to implment abstol=1e-9, reltol=1e-9
+    @time sol1 = solve(prob1, Rodas5P(); adaptive = false, dt = 1e-6) #TODO: fix BigFloat to implment abstol=1e-9, reltol=1e-9
     @time sol2 = solve(prob2, Rodas5P())
 
     # case 1: no negative pressure will only have gravity pulling mass back down