From 2d576cd57090076216437a8c33bff29ebe03c6d8 Mon Sep 17 00:00:00 2001 From: Venkateshprasad <32921645+ven-k@users.noreply.github.com> Date: Thu, 20 Mar 2025 06:20:31 +0000 Subject: [PATCH 1/2] docs: drop ref to non-existant `p_int` --- .../IsothermalCompressible/components.jl | 21 +++++++------------ .../IsothermalCompressible/sources.jl | 6 +++--- src/Hydraulic/IsothermalCompressible/utils.jl | 14 +++++-------- 3 files changed, 15 insertions(+), 26 deletions(-) diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index f7512bc48..669b1a713 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -1,12 +1,9 @@ """ - Cap(; p_int, name) + Cap(; name) Caps a hydraulic port to prevent mass flow in or out. -# Parameters: -- `p_int`: [Pa] initial pressure (set by `p_int` argument) - # Connectors: - `port`: hydraulic port """ @@ -14,7 +11,7 @@ Caps a hydraulic port to prevent mass flow in or out. @variables begin p(t), [guess = 0] - end + end @components begin port = HydraulicPort() @@ -28,12 +25,9 @@ Caps a hydraulic port to prevent mass flow in or out. end """ - Open(; p_int, name) + Open(; name) -Provides an "open" boundary condition for a hydraulic port such that mass flow `dm` is non-zero. This is opposite from an un-connected hydraulic port or the `Cap` boundary component which sets the mass flow `dm` to zero. - -# Parameters: -- `p_int`: [Pa] initial pressure (set by `p_int` argument) +Provides an "open" boundary condition for a hydraulic port such that mass flow `dm` is non-zero. This is opposite from an un-connected hydraulic port or the `Cap` boundary component which sets the mass flow `dm` to zero. # Connectors: - `port`: hydraulic port @@ -230,12 +224,11 @@ end @deprecate Pipe Tube """ - FlowDivider(;p_int, n, name) + FlowDivider(; 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: @@ -738,7 +731,7 @@ end """ SpoolValve2Way(reversible = false; p_s_int, p_a_int, p_b_int, p_r_int, m, g, x_int, Cd, d, name) -2-ways spool valve with 4 ports and spool mass. Fluid flow direction S → A and B → R when `x` is positive and S → B and A → R when `x` is negative. +2-ways spool valve with 4 ports and spool mass. Fluid flow direction S → A and B → R when `x` is positive and S → B and A → R when `x` is negative. # Parameters: - `p_s_int`: [Pa] initial pressure for `port_s` @@ -819,7 +812,7 @@ end Cd = 1e4, Cd_reverse = Cd, name) - + Actuator made of two DynamicVolumes connected in opposite direction with body mass attached. # Features: diff --git a/src/Hydraulic/IsothermalCompressible/sources.jl b/src/Hydraulic/IsothermalCompressible/sources.jl index 64b182f11..70939369c 100644 --- a/src/Hydraulic/IsothermalCompressible/sources.jl +++ b/src/Hydraulic/IsothermalCompressible/sources.jl @@ -1,12 +1,12 @@ """ - MassFlow(; name, p_int) + MassFlow(; name) Hydraulic mass flow input source # Connectors: - `port`: hydraulic port - - `dm`: real input + - `dm`: real input """ @mtkmodel MassFlow begin @@ -55,7 +55,7 @@ input pressure source # Connectors: - `port`: hydraulic port -- `p`: real input +- `p`: real input """ @mtkmodel Pressure begin diff --git a/src/Hydraulic/IsothermalCompressible/utils.jl b/src/Hydraulic/IsothermalCompressible/utils.jl index 425686d1f..85fbb247a 100644 --- a/src/Hydraulic/IsothermalCompressible/utils.jl +++ b/src/Hydraulic/IsothermalCompressible/utils.jl @@ -4,14 +4,10 @@ regPow(x, a, delta = 0.01) = x * (x * x + delta * delta)^((a - 1) / 2); regRoot(x, delta = 0.01) = regPow(x, 0.5, delta) """ - HydraulicPort(;p_int, name) + HydraulicPort(; name) Connector port for hydraulic components. -# Arguments: - -- `p_int`: [Pa] initial gauge pressure - # States: - `p`: [Pa] gauge total pressure - `dm`: [kg/s] mass flow @@ -38,7 +34,7 @@ end """ HydraulicFluid(; density = 997, bulk_modulus = 2.09e9, viscosity = 0.0010016, gas_density = 0.0073955, gas_pressure = -1000, n = 1, let_gas = 1, name) -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 (when `let_gas` is set to 1), this helps prevent pressures from going below the 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 (when `let_gas` is set to 1), this helps prevent pressures from going below the reference pressure. # Parameters: @@ -46,7 +42,7 @@ Fluid parameter setter for isothermal compressible fluid domain. Defaults given - `Β`: [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 -- `let_gas`: set to 1 to allow fluid to transition from liquid to gas (for density calculation only) +- `let_gas`: set to 1 to allow fluid to transition from liquid to gas (for density calculation only) - `ρ_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) """ @@ -80,12 +76,12 @@ f_turbulent(shape_factor, Re) = (shape_factor / 64) / (0.79 * log(Re) - 1.64)^2 """ friction_factor(dm, area, d_h, viscosity, shape_factor) -Calculates the friction factor ``f`` for fully developed flow in a tube such that ``Δp = f \\cdot \\rho \\frac{u^2}{2} \\frac{l}{d_h}`` where +Calculates the friction factor ``f`` for fully developed flow in a tube such that ``Δp = f \\cdot \\rho \\frac{u^2}{2} \\frac{l}{d_h}`` where - ``Δp``: [Pa] is the pressure difference over the tube length ``l`` - ``\\rho``: [kg/m^3] is the average fluid density - ``u``: [m/s] is the average fluid velocity -- ``l``: [m] is the tube length +- ``l``: [m] is the tube length The friction factor is calculated for laminar and turbulent flow with a transition region between Reynolds number 2000 to 3000. Turbulent flow equation is for smooth tubes, valid for the Reynolds number range up to 5e6. From 6a038b0a92dd41adc11e49828146116b351e92cc Mon Sep 17 00:00:00 2001 From: Venkateshprasad <32921645+ven-k@users.noreply.github.com> Date: Thu, 20 Mar 2025 06:21:44 +0000 Subject: [PATCH 2/2] style: fix format --- docs/src/tutorials/input_component.md | 30 ++++++++----------- src/Blocks/Blocks.jl | 3 +- src/Blocks/math.jl | 2 +- .../IsothermalCompressible/components.jl | 11 ------- .../IsothermalCompressible/sources.jl | 5 ---- src/Thermal/HeatTransfer/sources.jl | 6 ++-- src/Thermal/utils.jl | 3 +- test/Blocks/math.jl | 6 ++-- test/Blocks/test_analysis_points.jl | 15 ++++++---- test/Thermal/piston.jl | 5 ++-- test/Thermal/thermal.jl | 16 +++++----- 11 files changed, 44 insertions(+), 58 deletions(-) diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md index f48dd96b1..5d2843663 100644 --- a/docs/src/tutorials/input_component.md +++ b/docs/src/tutorials/input_component.md @@ -38,8 +38,7 @@ function MassSpringDamper(; name) @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0 @parameters m=10 k=1000 d=1 - eqs = [ - f ~ input.u + eqs = [f ~ input.u ddx * 10 ~ k * x + d * dx + f D(x) ~ dx D(dx) ~ ddx] @@ -52,10 +51,8 @@ function MassSpringDamperSystem(data, time; name) @named clk = ContinuousClock() @named model = MassSpringDamper() - eqs = [ - connect(src.input, clk.output) - connect(src.output, model.input) - ] + eqs = [connect(src.input, clk.output) + connect(src.output, model.input)] ODESystem(eqs, t, [], []; name, systems = [src, clk, model]) end @@ -100,11 +97,10 @@ using Plots function MassSpringDamper(; name) @named input = RealInput() - vars = @variables f(t) x(t)=0 dx(t) [guess=0] ddx(t) + vars = @variables f(t) x(t)=0 dx(t) [guess = 0] ddx(t) pars = @parameters m=10 k=1000 d=1 - eqs = [ - f ~ input.u + eqs = [f ~ input.u ddx * 10 ~ k * x + d * dx + f D(x) ~ dx D(dx) ~ ddx] @@ -117,10 +113,8 @@ function MassSpringDamperSystem(data, time; name) @named clk = ContinuousClock() @named model = MassSpringDamper() - eqs = [ - connect(model.input, src.output) - connect(src.input, clk.output) - ] + eqs = [connect(model.input, src.output) + connect(src.input, clk.output)] ODESystem(eqs, t; name, systems = [src, clk, model]) end @@ -143,6 +137,7 @@ plot(sol) ``` If we want to run a new data set, this requires only remaking the problem and solving again + ```@example parametrized_interpolation prob2 = remake(prob, p = [sys.src.data => ones(length(df.data))]) sol2 = solve(prob2) @@ -150,6 +145,7 @@ plot(sol2) ``` !!! note + Note that when changing the data, the length of the new data must be the same as the length of the original data. ## Custom Component with External Data @@ -224,7 +220,7 @@ using ModelingToolkitStandardLibrary.Blocks using OrdinaryDiffEq function System(; name) - src = SampledData(Float64, name=:src) + src = SampledData(Float64, name = :src) vars = @variables f(t)=0 x(t)=0 dx(t)=0 ddx(t)=0 pars = @parameters m=10 k=1000 d=1 @@ -238,7 +234,7 @@ function System(; name) end @named system = System() -sys = structural_simplify(system, split=false) +sys = structural_simplify(system, split = false) s = complete(system) dt = 4e-4 @@ -246,14 +242,14 @@ time = 0:dt:0.1 data1 = sin.(2 * pi * time * 100) data2 = cos.(2 * pi * time * 50) -prob = ODEProblem(sys, [], (0, time[end]); split=false, tofloat = false, use_union=true) +prob = ODEProblem(sys, [], (0, time[end]); split = false, tofloat = false, use_union = true) 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, build_initializeprob=false) + remake(prob; p, build_initializeprob = false) end prob1 = get_prob(data1) diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl index b4b18948c..55e1bd51a 100644 --- a/src/Blocks/Blocks.jl +++ b/src/Blocks/Blocks.jl @@ -10,7 +10,8 @@ using ModelingToolkit: getdefault, t_nounits as t, D_nounits as D export RealInput, RealInputArray, RealOutput, RealOutputArray, SISO include("utils.jl") -export Gain, Sum, MatrixGain, Feedback, Add, Add3, Product, Division, Power, Modulo, UnaryMinus, Floor, Ceil +export Gain, Sum, MatrixGain, Feedback, Add, Add3, Product, Division, Power, Modulo, + UnaryMinus, Floor, Ceil export Abs, Sign, Sqrt, Sin, Cos, Tan, Asin, Acos, Atan, Atan2, Sinh, Cosh, Tanh, Exp export Log, Log10 include("math.jl") diff --git a/src/Blocks/math.jl b/src/Blocks/math.jl index a787a9f43..d2ad76fc0 100644 --- a/src/Blocks/math.jl +++ b/src/Blocks/math.jl @@ -226,7 +226,7 @@ Output the exponential with base as the first input and exponent as second input output = RealOutput() end @equations begin - output.u ~ base.u ^ exponent.u + output.u ~ base.u^exponent.u end end diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl index 669b1a713..58d665963 100644 --- a/src/Hydraulic/IsothermalCompressible/components.jl +++ b/src/Hydraulic/IsothermalCompressible/components.jl @@ -8,7 +8,6 @@ Caps a hydraulic port to prevent mass flow in or out. - `port`: hydraulic port """ @mtkmodel Cap begin - @variables begin p(t), [guess = 0] end @@ -21,7 +20,6 @@ Caps a hydraulic port to prevent mass flow in or out. port.p ~ p port.dm ~ 0 end - end """ @@ -33,7 +31,6 @@ Provides an "open" boundary condition for a hydraulic port such that mass flow ` - `port`: hydraulic port """ @mtkmodel Open begin - @variables begin p(t), [guess = 0] dm(t), [guess = 0] @@ -47,7 +44,6 @@ Provides an "open" boundary condition for a hydraulic port such that mass flow ` port.p ~ p port.dm ~ dm end - end """ @@ -261,7 +257,6 @@ Reduces the flow from `port_a` to `port_b` by `n`. Useful for modeling parallel open.dm ~ dm_a - dm_b # extra flow dumps into an open port # port_b.dm ~ dm_b # divided flow goes to port_b end - end @component function ValveBase( @@ -355,7 +350,6 @@ Valve with `area` input and discharge coefficient `Cd` defined by https://en.wik end @mtkmodel VolumeBase begin - @structural_parameters begin Χ1 = 1 Χ2 = 1 @@ -385,7 +379,6 @@ end rho ~ full_density(port, port.p) port.dm ~ drho * vol * Χ1 + rho * area * dx * Χ2 end - end """ @@ -400,7 +393,6 @@ Fixed fluid volume. - `port`: hydraulic port """ @mtkmodel FixedVolume begin - @parameters begin vol end @@ -419,7 +411,6 @@ Fixed fluid volume. rho ~ full_density(port, port.p) port.dm ~ drho * vol end - end """ @@ -461,7 +452,6 @@ dm ────► │ │ area See also [`FixedVolume`](@ref), [`DynamicVolume`](@ref) """ @mtkmodel Volume begin - @structural_parameters begin direction = 1 end @@ -505,7 +495,6 @@ See also [`FixedVolume`](@ref), [`DynamicVolume`](@ref) @defaults begin rho => liquid_density(port) end - end """ diff --git a/src/Hydraulic/IsothermalCompressible/sources.jl b/src/Hydraulic/IsothermalCompressible/sources.jl index 70939369c..cb032a109 100644 --- a/src/Hydraulic/IsothermalCompressible/sources.jl +++ b/src/Hydraulic/IsothermalCompressible/sources.jl @@ -9,7 +9,6 @@ Hydraulic mass flow input source - `dm`: real input """ @mtkmodel MassFlow begin - @components begin port = HydraulicPort() dm = RealInput() @@ -32,7 +31,6 @@ Fixed pressure source - `port`: hydraulic port """ @mtkmodel FixedPressure begin - @parameters begin p end @@ -44,7 +42,6 @@ Fixed pressure source @equations begin port.p ~ p end - end @deprecate Source FixedPressure @@ -58,7 +55,6 @@ input pressure source - `p`: real input """ @mtkmodel Pressure begin - @components begin port = HydraulicPort() p = RealInput() @@ -67,6 +63,5 @@ input pressure source @equations begin port.p ~ p.u end - end @deprecate InputSource Pressure diff --git a/src/Thermal/HeatTransfer/sources.jl b/src/Thermal/HeatTransfer/sources.jl index 1057f34d8..cdd279479 100644 --- a/src/Thermal/HeatTransfer/sources.jl +++ b/src/Thermal/HeatTransfer/sources.jl @@ -28,9 +28,9 @@ the component FixedHeatFlow is connected, if parameter `Q_flow` is positive. end @equations begin - port.Q_flow ~ ifelse(alpha == 0.0, - -Q_flow, # Simplified equation when alpha is 0 - -Q_flow * (1 + alpha * (port.T - T_ref))) + port.Q_flow ~ ifelse(alpha == 0.0, + -Q_flow, # Simplified equation when alpha is 0 + -Q_flow * (1 + alpha * (port.T - T_ref))) end end diff --git a/src/Thermal/utils.jl b/src/Thermal/utils.jl index 1aa5c53bf..22c379a07 100644 --- a/src/Thermal/utils.jl +++ b/src/Thermal/utils.jl @@ -1,4 +1,5 @@ -@connector function HeatPort(; T = nothing, T_guess = 273.15 + 20, Q_flow = nothing, Q_flow_guess = 0.0, name) +@connector function HeatPort(; + T = nothing, T_guess = 273.15 + 20, Q_flow = nothing, Q_flow_guess = 0.0, name) pars = @parameters begin T_guess = T_guess Q_flow_guess = Q_flow_guess diff --git a/test/Blocks/math.jl b/test/Blocks/math.jl index 5ee902604..262a99236 100644 --- a/test/Blocks/math.jl +++ b/test/Blocks/math.jl @@ -180,7 +180,7 @@ end @testset "Modulo" begin @named c1 = Ramp(height = 2, duration = 1, offset = 1, start_time = 0, smooth = false) - @named c2 = Constant(; k = 1) + @named c2 = Constant(; k = 1) @named modl = Modulo(;) @named model = ODESystem( [ @@ -194,7 +194,7 @@ end sol = solve(prob, Rodas4()) @test isequal(unbound_inputs(sys), []) @test sol.retcode == Success - @test sol[modl.remainder.u] ≈ mod.(2 * sol.t,1) + @test sol[modl.remainder.u] ≈ mod.(2 * sol.t, 1) end @testset "UnaryMinus" begin @@ -213,7 +213,7 @@ end sol = solve(prob, Rodas4()) @test isequal(unbound_inputs(sys), []) @test sol.retcode == Success - @test sol[minu.output.u] ≈ - sin.(2 * pi * sol.t) + @test sol[minu.output.u] ≈ -sin.(2 * pi * sol.t) end @testset "Floor" begin diff --git a/test/Blocks/test_analysis_points.jl b/test/Blocks/test_analysis_points.jl index c92e34e24..4d7d2bfb8 100644 --- a/test/Blocks/test_analysis_points.jl +++ b/test/Blocks/test_analysis_points.jl @@ -266,7 +266,8 @@ Sinner = sminreal(ss(get_sensitivity(sys_inner, :u)[1]...)) Souter = sminreal(ss(get_sensitivity(sys_outer, sys_outer.sys_inner.u)[1]...)) -Sinner2 = sminreal(ss(get_sensitivity(sys_outer, sys_outer.sys_inner.u, loop_openings = [:y2])[1]...)) +Sinner2 = sminreal(ss(get_sensitivity( + sys_outer, sys_outer.sys_inner.u, loop_openings = [:y2])[1]...)) @test Sinner.nx == 1 @test Sinner == Sinner2 @@ -332,7 +333,8 @@ eqs = [connect(r.output, F.input) connect(F.output, sys_inner.add.input1)] sys_outer = ODESystem(eqs, t, systems = [F, sys_inner, r], name = :outer) -matrices, _ = get_sensitivity(sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output]) +matrices, _ = get_sensitivity( + sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output]) Ps = tf(1, [1, 1]) |> ss Cs = tf(1) |> ss @@ -346,7 +348,8 @@ So = CS.feedback(1, Ps * Cs) @test tf(G[1, 2]) ≈ tf(-CS.feedback(Cs, Ps)) @test tf(G[2, 1]) ≈ tf(CS.feedback(Ps, Cs)) -matrices, _ = get_comp_sensitivity(sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output]) +matrices, _ = get_comp_sensitivity( + sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output]) G = CS.ss(matrices...) |> sminreal Ti = CS.feedback(Cs * Ps) @@ -369,13 +372,15 @@ L = CS.ss(matrices...) |> sminreal @test tf(L[1, 1]) ≈ -tf(Ps * Cs) # Calling looptransfer like below is not the intended way, but we can work out what it should return if we did so it remains a valid test -matrices, _ = get_looptransfer(sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output]) +matrices, _ = get_looptransfer( + sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output]) L = CS.ss(matrices...) |> sminreal @test tf(L[1, 1]) ≈ tf(0) @test tf(L[2, 2]) ≈ tf(0) @test sminreal(L[1, 2]) ≈ ss(-1) @test tf(L[2, 1]) ≈ tf(Ps) -matrices, _ = linearize(sys_outer, [sys_outer.inner.plant_input], [sys_outer.inner.plant_output]) +matrices, _ = linearize( + sys_outer, [sys_outer.inner.plant_input], [sys_outer.inner.plant_output]) G = CS.ss(matrices...) |> sminreal @test tf(G) ≈ tf(CS.feedback(Ps, Cs)) diff --git a/test/Thermal/piston.jl b/test/Thermal/piston.jl index c8b2c2563..9b915f719 100644 --- a/test/Thermal/piston.jl +++ b/test/Thermal/piston.jl @@ -42,6 +42,7 @@ using ModelingToolkitStandardLibrary.Blocks # The initial value doesn't add up to absolute zero, while the rest do. To avoid # tolerance on the latter, the test is split in two parts. @test sol[piston.gas.Q_flow][1] + sol[piston.coolant.Q_flow][1]≈0 atol=1e-6 - @test sol[piston.gas.Q_flow][2:end] + sol[piston.coolant.Q_flow][2:end] ≈ - zeros(length(sol) - 1) atol = 1e-6 + @test sol[piston.gas.Q_flow][2:end] + + sol[piston.coolant.Q_flow][2:end]≈ + zeros(length(sol) - 1) atol=1e-6 end diff --git a/test/Thermal/thermal.jl b/test/Thermal/thermal.jl index 9513aa260..43137d9d5 100644 --- a/test/Thermal/thermal.jl +++ b/test/Thermal/thermal.jl @@ -158,20 +158,18 @@ end @testset "FixedHeatFlow with alpha=0.0 test" begin @mtkmodel TestModel begin - @components begin - temp = FixedTemperature(T=300) - heatflow = FixedHeatFlow(Q_flow=-1.0) - wall = ThermalResistor(R=1) + temp = FixedTemperature(T = 300) + heatflow = FixedHeatFlow(Q_flow = -1.0) + wall = ThermalResistor(R = 1) end - + @equations begin connect(temp.port, wall.port_a) connect(wall.port_b, heatflow.port) end - end - + @info "Building a FixedHeatFlow with alpha=0.0" @mtkbuild test_model = TestModel() prob = ODEProblem(test_model, Pair[], (0, 10.0)) @@ -180,5 +178,5 @@ end heat_flow = sol[test_model.heatflow.port.Q_flow] @test SciMLBase.successful_retcode(sol) # Ensure the simulation is successful - @test all(isapprox.(heat_flow, 1.0, rtol=1e-6)) # Heat flow value should be equal to the fixed value defined -end \ No newline at end of file + @test all(isapprox.(heat_flow, 1.0, rtol = 1e-6)) # Heat flow value should be equal to the fixed value defined +end