Skip to content

Commit

Permalink
Adding acausal equations
Browse files Browse the repository at this point in the history
  • Loading branch information
wallytutor committed Oct 17, 2024
1 parent 92d043c commit 19b51bd
Show file tree
Hide file tree
Showing 3 changed files with 141 additions and 114 deletions.
123 changes: 9 additions & 114 deletions docs/src/Notebooks/Pluto/darcy-weisbach.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@ begin
using ModelingToolkit
using NonlinearSolve
using PlutoUI: TableOfContents
using Symbolics

TableOfContents()
end
Expand Down Expand Up @@ -256,112 +255,6 @@ md"""
## Experimental
"""

# ╔═╡ db5fe5d7-e033-43ef-92d3-1dcbf3d64c17
begin
tuplefy(d) = NamedTuple{Tuple(keys(d))}(values(d))

########

struct UnsolvableAlgebraicModel <: Exception
msg::String

function UnsolvableAlgebraicModel(vars, pars)
return new(
"Unable to solve algebraic model containing [$(join(vars, ", "))] " *
"by providing only [$(join(pars, ", "))]! Please, revise the " *
"setup of the system!"
)
end
end

function Base.showerror(io::IO, err::UnsolvableAlgebraicModel)
print(io, "UnsolvableAlgebraicModel: ")
print(io, err.msg)
end

########

struct AlgebraicModel
vars::NamedTuple
eqns::Vector{Equation}
pars::Dict{Symbol, Float64}

function AlgebraicModel(vars, eqns, pars)
vars = Dict(map(v->(Symbol(v), v), vars))
return new(tuplefy(vars), eqns, pars)
end
end

function replacement_solve(obj::AlgebraicModel)
subs = map(v->(v, get(obj.pars, Symbol(v), v)), values(obj.vars))
eqns = substitute(obj.eqns, Dict(subs))
return eqns
end

function structural_solve(obj::AlgebraicModel)
eqns = replacement_solve(obj)

has_pars = length(obj.pars) > 0

if has_pars && length(obj.vars) - length(obj.pars) != length(eqns)
throw(UnsolvableAlgebraicModel(obj.vars, obj.pars))
end

@mtkbuild ns = NonlinearSystem(eqns, [values(obj.vars)...], [])
prob = NonlinearProblem(ns, [], [])
sol = solve(prob)

solved_for = map(e->e.lhs, observed(ns))

pars = tuplefy(obj.pars)
vars = NamedTuple(map(n->(Symbol(n), sol[n]), solved_for))
return merge(pars, vars)
end

########

function perfect_gas_constants(; pars...)
vars = @variables begin
k, [bounds = (1, Inf)]
R, [bounds = (0, Inf)]
cp, [bounds = (0, Inf)]
cv, [bounds = (0, Inf)]
end

eqns = [
0 ~ cp - k * cv,
0 ~ cv - R / (k - 1)
]

return AlgebraicModel(vars, eqns, pars)
end

function perfect_gas(; pars...)
model = perfect_gas_constants(;)

vars = @variables begin
p, [bounds = (0, Inf)]
ρ, [bounds = (0, Inf)]
T, [bounds = (0, Inf)]
h, [bounds = (0, Inf)]
s, [bounds = (0, Inf)]
end

R = model.vars.R
cp = model.vars.cp

vars = [vars..., values(model.vars)...]
eqns = [
0 ~ p - ρ * R * T,
0 ~ h - cp * T,
0 ~ s - cp * log(T) + R * log(p),
model.eqns...
]

return AlgebraicModel(vars, eqns, pars)
end
end

# ╔═╡ be031284-d872-4e1c-b43d-53038b21f8e4
let
model = perfect_gas_constants(; R=208, k=1.67)
Expand All @@ -373,8 +266,8 @@ let
R = 208.0
k = 1.67

model1 = perfect_gas(; R, k, p=1.70e+06, ρ=18.0)
model2 = perfect_gas(; R, k, p=2.48e+05, T=400.0)
model1 = perfect_gas_model(; R, k, p=1.70e+06, ρ=18.0)
model2 = perfect_gas_model(; R, k, p=2.48e+05, T=400.0)

sol1 = structural_solve(model1)
sol2 = structural_solve(model2)
Expand All @@ -389,9 +282,12 @@ end
let
R = 287.0
k = 1.4

p1 = P_REF + 3550.0
p2 = P_REF

model1 = perfect_gas(; R, k, p=P_REF+3550.0, T=300.0)
model2 = perfect_gas(; R, k, p=P_REF, T=300.0)
model1 = perfect_gas_model(; R, k, p=p1, T=300.0)
model2 = perfect_gas_model(; R, k, p=p2, T=300.0)

sol1 = structural_solve(model1)
sol2 = structural_solve(model2)
Expand Down Expand Up @@ -420,8 +316,8 @@ end
# ╟─360e9495-3384-477b-a211-ffa8bb3d1fac
# ╟─00812f04-a14c-434f-a0ac-bf0943976a0c
# ╟─795203f6-f46b-4a47-8cd4-f50b2dc7ef38
# ╠═882d83ed-65e5-4b24-b4d2-305e5f63baf3
# ╠═60985445-5a27-41fd-9956-c88f47704274
# ╟─882d83ed-65e5-4b24-b4d2-305e5f63baf3
# ╟─60985445-5a27-41fd-9956-c88f47704274
# ╟─a02cd203-8122-444d-a0cc-af32ce1ffd5c
# ╟─bd79f2ef-db5e-442b-98eb-a209820cbefe
# ╠═1f115f6f-df6d-4134-8269-167736e3c238
Expand All @@ -431,7 +327,6 @@ end
# ╟─3357de42-e48c-4e4d-a7df-58d916a1c9b5
# ╟─8fcaed37-d179-449c-8d43-8d0ad9a6c35e
# ╟─6a6ad275-90d9-47d0-8b5e-a3d20efa9b0a
# ╠═db5fe5d7-e033-43ef-92d3-1dcbf3d64c17
# ╠═be031284-d872-4e1c-b43d-53038b21f8e4
# ╠═80dd3540-f5fb-4c2f-b6d5-b39ad599b8d5
# ╠═3e85af5c-7160-4b60-9369-dfa9dad34f7e
4 changes: 4 additions & 0 deletions src/wallytoolbox/thermochemistry.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
import Base: +
import Base: *
import Base: String
import Base: showerror

import DataFrames
import YAML
Expand All @@ -13,11 +14,13 @@ using CairoMakie
using Distributions: Weibull
using Distributions: cdf, scale, mean, params
using DocStringExtensions: TYPEDFIELDS
using ModelingToolkit
using Polynomials: AbstractPolynomial
using Printf
using Polynomials: Polynomial, LaurentPolynomial, integrate
using Roots
using SteamTables: SpecificV
using Symbolics

##############################################################################
# CONSTANT AND CONFIGURATION
Expand All @@ -44,6 +47,7 @@ end
include("thermochemistry/chemistry.jl")
include("thermochemistry/combustion.jl")
include("thermochemistry/thermodata.jl")
include("thermochemistry/acausal.jl")

##############################################################################
# EOF
Expand Down
128 changes: 128 additions & 0 deletions src/wallytoolbox/thermochemistry/acausal.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,128 @@
##############################################################################
# ACAUSAL
##############################################################################

export AlgebraicModel
export structural_solve
export perfect_gas_constants
export perfect_gas_model

##############################################################################
# EXCEPTIONS
##############################################################################

struct UnsolvableAlgebraicModel <: Exception
msg::String

function UnsolvableAlgebraicModel(vars, pars)
svars, spars = join(vars, ", "), join(pars, ", ")

return new(
"Unable to solve algebraic model containing [$(svars)] " *
"by providing only [$(spars)]! Please, revise the " *
"setup of the system!"
)
end
end

function Base.showerror(io::IO, err::UnsolvableAlgebraicModel)
print(io, "UnsolvableAlgebraicModel: ")
print(io, err.msg)
end

##############################################################################
# ALGEBRAIC MODEL
##############################################################################

struct AlgebraicModel
vars::NamedTuple
eqns::Vector{Equation}
pars::Dict{Symbol, Float64}

function AlgebraicModel(vars, eqns, pars)
vars = Dict(map(v->(Symbol(v), v), vars))
return new(tuplefy(vars), eqns, pars)
end
end

function structural_solve(obj::AlgebraicModel)
eqns = replacement_solve(obj)

if has_pars(obj) && has_dofs(obj, eqns)
throw(UnsolvableAlgebraicModel(obj.vars, obj.pars))
end

@mtkbuild ns = NonlinearSystem(eqns, [values(obj.vars)...], [])
sol = solve(NonlinearProblem(ns, [], []))

solved_for = map(e->e.lhs, observed(ns))

pars = tuplefy(obj.pars)
vars = NamedTuple(map(n->(Symbol(n), sol[n]), solved_for))
return merge(pars, vars)
end

function replacement_solve(obj::AlgebraicModel)
# Retrieve names provided in object `pars` (parameters).
subs = map(v->(v, get(obj.pars, Symbol(v), v)), values(obj.vars))

# Replace at once all identified symbols for solving.
return substitute(obj.eqns, Dict(subs))
end

function has_pars(obj::AlgebraicModel)
return length(obj.pars) > 0
end

function has_dofs(obj::AlgebraicModel, eqns)
return length(obj.vars) - length(obj.pars) != length(eqns)
end

##############################################################################
# APPLICATIONS
##############################################################################

function perfect_gas_constants(; pars...)
vars = @variables begin
k, [bounds = (1, Inf)]
R, [bounds = (0, Inf)]
cp, [bounds = (0, Inf)]
cv, [bounds = (0, Inf)]
end

eqns = [
0 ~ cp - k * cv,
0 ~ cv - R / (k - 1)
]

return AlgebraicModel(vars, eqns, pars)
end

function perfect_gas_model(; pars...)
model = perfect_gas_constants(;)

vars = @variables begin
p, [bounds = (0, Inf)]
ρ, [bounds = (0, Inf)]
T, [bounds = (0, Inf)]
h, [bounds = (0, Inf)]
s, [bounds = (0, Inf)]
end

R = model.vars.R
cp = model.vars.cp

vars = [vars..., values(model.vars)...]
eqns = [
0 ~ p - ρ * R * T,
0 ~ h - cp * T,
0 ~ s - cp * log(T) + R * log(p),
model.eqns...
]

return AlgebraicModel(vars, eqns, pars)
end

##############################################################################
# EOF
##############################################################################

0 comments on commit 19b51bd

Please sign in to comment.