Fuse Ready

.gitignore

@@ -1 +1,7 @@
 /Manifest.toml
+src/Package_testing.jl
+src/Playground/Gtesting.jl
+src/Playground/Package_testing.jl
+src/Playground/Package_testing2.jl
+src/Playground/pcktst.jl
+src/Playground/Workflows.jl

Project.toml

@@ -5,23 +5,38 @@ version = "0.0.1"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 CoolProp = "e084ae63-2819-5025-826e-f8e611a84251"
 Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
 GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231"
 GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
+Ipopt_jll = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+LayeredLayouts = "f4a74d36-062a-4d48-97cd-1356bad1de4e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
+NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
 NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
+OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -31,6 +46,7 @@ PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

src/DeadTime/02-nodes.jl

@@ -284,16 +284,16 @@ liq_properties(n::component_nodes) = liq_properties(n.working_fluid;P = n.P, T =
 gas_properties(g::component_nodes) = gas_properties(g.working_fluid;P=g.P,T=g.T)
 
 
-function gastype(working_fluid::Symbol)
+# function gastype(working_fluid::Symbol)
 
-    if working_fluid in [:air, :helium] 
-        return :gas
-    else
-        return :liq
-    end
+#     if working_fluid in [:air, :helium] 
+#         return :gas
+#     else
+#         return :liq
+#     end
 
 
-end
+# end
 
 function calculate_fluid_properties(n::component_nodes)
     fluid =  n.working_fluid

src/ODE_Systems/01-Incompressible.jl

@@ -4,13 +4,18 @@ using NonlinearSolve, Printf
 using Plots
 @variables t
 Logging.disable_logging(Logging.Warn)
-include("03-MTK_UTILS.jl")
+# include("03-MTK_UTILS.jl")
 # data from [2]
 cppb(x)     = (0.195 - 9.116e-6 .* x) .* 1000  # J/kgK
 vpb(x)      = 1 / (10520.35 - 1.19051 .* x)      # m³/kg
 spb(x)      = 0.195*1e3*log(x) - 9.11e-6*1000*x - 1210.316002277804
 hpb(x,y)      = cppb(x)*x +  vpb(x) * y * 1e5 
 
+Lcpfunc(x) = cppb(x)
+Lvfunc(x) = vpb(x)
+Lsfunc(x) = spb(x)
+Lhfunc(x,y) = hpb(x,y)
+
 @variables x
 @register_symbolic Lcpfunc(x)
 @register_symbolic Lvfunc(x)
@@ -22,16 +27,17 @@ propDict = Dict(Lcpfunc => cppb,
             Lhfunc => hpb)
 
 
-@connector function IncompressiblePin(; name, Pdef = 50, Tdef = 555)
-    across_var  = @variables  P(t)=Pdef T(t)=Tdef s(t)=0.0 cp(t)=187 v(t)=.001
-    thru_var    = @variables  ṁ(t)=0.0 Φ(t)=0.0                     # mass flow and energy flow
+@connector function IncompressiblePin(; name, Pdef = 50, Tdef = 555, ṁdef = 0.0)
+    across_var  = @variables  P(t)=Pdef T(t)=Tdef s(t)=1.0 cp(t)=187 v(t)=.001
+    thru_var    = @variables  ṁ(t)=ṁdef Φ(t)=1.0                     # mass flow and energy flow
 
     eq = [
         cp ~ Lcpfunc(T)
         v ~ Lvfunc(T)
         s ~ Lsfunc(T)
     ]
-    ODESystem(eq, t, [across_var...,thru_var...], []; name = name, defaults = [P => Pdef, T => Tdef, ṁ =>1.0, Φ => 0.0])
+    sts = [T, P, ṁ, cp,s,v, Φ]
+    ODESystem(eq, t, sts, []; name = name, defaults = [P => Pdef,T=>Tdef, ṁ =>ṁdef, Φ => 0.0])
 end
 
 @connector function WorkPin(; name)
@@ -79,6 +85,20 @@ end
         extend(ODESystem(eqs,t,[],ps; name = name), oneport)
 end
 
+@component function SetPressure2(;name, P = 0.1)
+    @named p = IncompressiblePin(Pdef = P)
+    @named n = IncompressiblePin(Pdef = P)
+    ps = @parameters P = P 
+    eqs = [
+        p.P ~ P
+        n.P ~ p.P
+        n.T ~ p.T
+        0 ~ p.Φ + n.Φ             # conservation of energy
+        0 ~ p.ṁ + n.ṁ
+    ]
+    ODESystem(eqs, t,[], ps; name = name, systems = [p,n])
+end
+
 @component function SetTemperature(;name, T = 300)
     @named oneport = IncompressibleOnePort()
     @unpack ΔT, ΔP, LHS, p = oneport
@@ -135,12 +155,30 @@ end
     extend(ODESystem(eqs,t,[],ps;systems = [w], name = name), oneport)
 end
 
-@component function IncompressibleHeaatTransfer(; name)
+@component function PassiveIncompressiblePump2(; name, η = 1.0)
+    @named p = IncompressiblePin()
+    @named n = IncompressiblePin()
+    @named w = WorkPin()
+
+    ps = @parameters η = η
+
+    eqs = [
+        0 ~ p.ṁ + n.ṁ #p.ṁ ~ n.ṁ                               # conservation of mass
+        w.Ẇ ~ p.ṁ * p.v * (n.P - p.P) * 1e5 / η
+        w.Ẇ ~ p.Φ + n.Φ
+        n.T  ~ p.T + (p.v * (n.P - p.P) * 1e5 / η - p.v*(n.P - p.P)*1e5)/p.cp
+    ]
+    ODESystem(eqs,t,[],ps; name = name, systems = [p,n,w])
+end
+
+
+
+@component function IncompressibleHeatTransfer(; name)
     @named p = IncompressiblePin()
     @named n = IncompressiblePin()
     @named q = HeatTransferPin()
 
-    st = @variables Q̇(t)=0.0 C(t)=5192
+    st = @variables Q̇(t)=0.0 C(t)=187
 
     eqs = [
         0 ~ p.ṁ + n.ṁ           # p.ṁ ~ n.ṁ                               # conservation of mass
@@ -153,15 +191,35 @@ end
     ODESystem(eqs,t,[Q̇,C],[]; name = name, systems = [p,n,q], defaults = [Q̇ => 0.0, C => 187])
 end
 
+@component function FlowControlIncompressibleHeatTransfer(; name, ΔP = 0.0, Tout = 1000.0)
+    @named p = IncompressiblePin()
+    @named n = IncompressiblePin()
+    @named q = HeatTransferPin()
+
+    st = @variables Q̇(t)=0.0 C(t)=187 ṁ(t)=1.0
+    ps = @parameters ΔP = ΔP Tout = Tout
+    eqs = [
+        n.T ~ Tout
+        ṁ ~ q.Q̇ /((n.T - p.T)*p.cp)           # q = mcp*ΔT
+        p.ṁ ~ ṁ
+        0 ~ p.ṁ + n.ṁ           # p.ṁ ~ n.ṁ                               # conservation of mass
+        q.Q̇ ~ p.Φ + n.Φ         # conservation of energy            
+        C ~ ṁ * p.cp          # duty
+        0 ~ q.Q̇ - Q̇ 
+        n.P ~ p.P - ΔP
+    ]
+    ODESystem(eqs,t,[Q̇,C,ṁ],ps; name = name, systems = [p,n,q], defaults = [Q̇ => 0.0, C => 187])
+end
+
 @component function SinglePortReservoir(;name, P = 0.1, T = 300)
     @named n = IncompressiblePin(Pdef = P, Tdef = T)
-    ps = @parameters P=P 
+    ps = @parameters P=P T=T
     eqs = [
         n.P ~ P
         n.T ~ T
         n.Φ ~ 0
     ]
-    ODESystem(eqs, t,[], ps; name = name, systems = [n], defaults = [P => 0.1, T =>300])
+    ODESystem(eqs, t,[], ps; name = name, systems = [n])
 end
 
 @component function TwoPortReservoir(;name, P = 0.1, T = 300)
@@ -204,6 +262,18 @@ end
     ODESystem(eqs, t, sts, []; name = name, systems = [p,n,q], defaults = [Q̇ => 0])
 end
 
+@component function ReliefElement(;name)
+    @named p = IncompressiblePin()
+    @named n = IncompressiblePin()
+    @named q = HeatTransferPin()
+    # 0 variables, 3 equations for pin
+    eqs = [
+        n.ṁ + p.ṁ ~ 0           # has to be negative
+        q.Q̇ ~ p.Φ + n.Φ         # conservation of energy
+        ]
+    ODESystem(eqs, t,[], []; name = name, systems = [p,n,q])
+end
+
 """
     incompressible_connect(pins...)
     Adds mass and energy balance eqs to all nodes
@@ -221,6 +291,7 @@ function incompressible_connect(pins...)
     end
     return eqs
 end
+
 function showsol(c,sol)
     for cel in c
         cvec = []

src/ODE_Systems/01-MultiPhase.jl

@@ -4,7 +4,7 @@ using Graphs, Plots, GraphRecipes
 using DifferentialEquations
 @variables t
 Logging.disable_logging(Logging.Warn)
-include("03-MTK_UTILS.jl")
+# include("03-MTK_UTILS.jl")
 ##
 #       Need todo: Coolprops, validation
 #               Connect work nodes - master energy balance?
@@ -246,7 +246,7 @@ end
         n.h  ~ stm_hsatfunc(P)
         n.Φ ~ 0
     ]
-    ODESystem(eqs, t,[], ps; name = name, systems = [n], defaults = [P => 0.1])
+    ODESystem(eqs, t,[], ps; name = name, systems = [n])
 end
 
 @component function MultiPhaseGnd(;name, P = 0.1)
@@ -256,7 +256,7 @@ end
         n.P ~ P
         n.Φ ~ 0
     ]
-    ODESystem(eqs, t,[], ps; name = name, systems = [n], defaults = [P => 0.1])
+    ODESystem(eqs, t,[], ps; name = name, systems = [n])
 end
 
 @component function SuperHeatedReservoir(;name, P = 150, T = 600)
@@ -267,7 +267,7 @@ end
         0 ~ stm_hptfunc(P,T) + n.h
         n.Φ ~ -n.ṁ * n.h
     ]
-    ODESystem(eqs, t,[], ps; name = name, systems = [n], defaults = [P => 0.1, T => T])
+    ODESystem(eqs, t,[], ps; name = name, systems = [n])
 end
 
 @component function ioReservoir(;name, P = 0.1, fixboth = false)
@@ -284,7 +284,7 @@ end
         push!(eqs,p.P ~ P)
         # push!(eqs,p.h ~ n.h)
     end
-    ODESystem(eqs, t,[], ps; name = name, systems = [n,p], defaults = [P => 0.1 ])
+    ODESystem(eqs, t,[], ps; name = name, systems = [n,p])
 end
 
 @component function TwoPortReservoir(;name, P = 0.1)
@@ -300,7 +300,7 @@ end
         0 ~ n.ṁ + p.ṁ                                                           # 1/density = specific volume
         0 ~ n.Φ + p.Φ
     ]
-    ODESystem(eqs, t,[], ps; name = name, systems = [n,p], defaults = [P => 0.1 ])
+    ODESystem(eqs, t,[], ps; name = name, systems = [n,p])
 end
 
 @component function ContinuityReservoir2(;name)
@@ -394,14 +394,13 @@ end
     @named p = BasicSteamPin()
     @named n = BasicSteamPin()
     @named w = WorkPin()
-
     ps = @parameters η = η P=Pout
 
     eqs = Equation[
         w.Ẇ ~ p.Φ + n.Φ                         # conservation of energy
         0 ~ p.ṁ + n.ṁ
         n.h  ~ p.h + p.v *1e5* (n.P - p.P) / η  # work, multiply by 100 to get to kPa
-        0  ~ p.ṁ * (n.h - p.h)  + w.Ẇ 
+        w.Ẇ  ~ p.ṁ * (n.h - p.h)
     ]
 
     if setpressure
@@ -412,7 +411,7 @@ end
         eqs = vcat(eqs,p.h  ~ stm_hsatfunc(p.P))
     end
     # extenda([ODESystem(eqs, t,[], ps; name = name, defaults = [η => 0.6, P => Pout]), p,n,w])
-    ODESystem(eqs, t,[], ps; name = name, systems = [p,n,w], defaults = [η => .6, P => 10])
+    ODESystem(eqs, t,[], ps; name = name, systems = [p,n,w])
 end
 
 @component function Splitter(;name)
@@ -447,7 +446,7 @@ end
         w.Ẇ ~ Ẇ
         0 ~ p.ṁ + n.ṁ      
         n.h  ~ p.h - (p.h-stm_hpsfunc(n.P,p.s))*η
-        0 ~ p.ṁ * (n.h - p.h) + w.Ẇ
+        w.Ẇ ~ p.ṁ * (n.h - p.h)
     ]
 
     if setpressure
@@ -485,7 +484,7 @@ end
 
     eqs = vcat(split_connect,hp_connect,lp_connect)
 
-    compose(ODESystem(eqs, t,[], ps; name = name, systems = [yn,zn,p], defaults = [η => 1.0, Py => 10, Pz => 0.1]), hp,lp)
+    compose(ODESystem(eqs, t,[], ps; name = name, systems = [yn,zn,p]), hp,lp)
     # ODESystem(eqs, t,[], ps; name = name, systems = [p,hp,lp], defaults = [η => 1.0, Py => 10, Pz => 0.1])
     # extend(ODESystem(eqs, t,sts, ps; name = name, systems = [hp,lp], defaults = [η => 1.0 Py => 10 Pz => 0.1]),split)
 end
@@ -503,7 +502,7 @@ end
         n.h ~ p.h + q.Q̇/(p.ṁ)
         n.P ~ p.P
     ]
-    ODESystem(eqs,t,[Q̇,C],[]; name = name, systems = [p,n,q], defaults = [Q̇ => 0.0, C => 400])
+    ODESystem(eqs,t,[Q̇,C],[]; name = name, systems = [p,n,q])
 end
 
 @component function TunableSteamHeatTransfer(;name, Q̇in = 150e6)
@@ -576,6 +575,21 @@ end
     ODESystem(eqs, t,sts, []; name = name, systems = [p,n,q], defaults = [Q̇ => 0.0])
 end
 
+@component function ReliefElement(;name, pressurecontrol = false)
+    @named p = BasicSteamPin()
+    @named n = BasicSteamPin()
+    @named q = HeatTransferPin()
+    # 0 variables, 3 equations for pin
+    eqs = [
+        n.ṁ + p.ṁ ~ 0           # has to be negative
+        q.Q̇ ~ p.Φ + n.Φ         # conservation of energy
+        ]
+    if pressurecontrol
+        eqs = vcat(eqs, n.P ~ p.P)
+    end
+    ODESystem(eqs, t,[], []; name = name, systems = [p,n,q])
+end
+
 @component function OpenFeedwaterHeater(;name)
     # flows x and y are the inlets
     @named p1 = BasicSteamPin()
@@ -596,7 +610,7 @@ end
         0 ~ n.ṁ + p1.ṁ + p2.ṁ
         0 ~ n.Φ + (p1.Φ + p2.Φ)
     ]
-    ODESystem(eqs, t,sts, []; name = name, systems = [p1,p2,n], defaults = [yfrac => 0.5], continuous_events)
+    ODESystem(eqs, t,sts, []; name = name, systems = [p1,p2,n], continuous_events)
 end
 
 @component function MixingChamber(;name)