ModelingToolkit KeyError

Specific Domains Modelling & Simulations
, ,
1

I’m trying to model powder burning in small chamber by the following code.
The problem is that I need to stop integration when there is more gas then chamber volume. I try to do it by adding callback to bal_eqs. But I get KeyError: key bal_eqs₊Wg(t) not found then I try to create ODEProblem. May be I need to add this terminate condition to connected ODESystem.Can anybody help me or should I write a GitHub Issue?

using ModelingToolkit,Plots, DifferentialEquations
function terminate_affect!(integ, u, p, ctx)
    terminate!(integ)
end
begin
	@variables t
	D = Differential(t)
end
function powder_burning(ν=0.75, p_s=50e6;name)
	@parameters κ1 λ1 μ1 κ2 λ2 μ2 Jk zk
	@variables t ψ(t) z(t) pressure(t)

	dψ = D(ψ) ~ ifelse(
			z < 1.0,
			κ1 * (1. + 2λ1 * z + 3μ1 * z^2) * D(z),
			κ2 * (1. + 2λ2 * (z-1) + 3μ2 * (z-1)^2) * D(z)
		)

	dz = ifelse(
		pressure > p_s,
		p_s^(1-ν) * pressure^ν / Jk,
		pressure / Jk
	)

	dz = D(z) ~ ifelse(
			z < zk,
			dz,
			0
		)
	ODESystem([dψ, dz]; name)
end
function powder_contribution(;name)
	@parameters f, θ, α, T_b, ω, ρ
	@variables t, dE(t), dW(t), θ_num(t), θ_denum(t), ψ(t)

	eqs = [
		dE ~ f * ω * ψ
		dW ~ (ω * (1 - ψ) / ρ) - ψ * α
		θ_num ~ f * ω * ψ / T_b
		θ_denum ~ (f * ω * ψ) / (T_b * θ)
		]
	ODESystem(eqs; name)
end
begin
	@named p1_b = powder_burning()
	@named p2_b = powder_burning()
end
begin
	@named p1_c = powder_contribution()
	@named p2_c = powder_contribution()
end
begin
	@variables p(t) E(t) Wg(t)
	@parameters W0
end
eqs = [
	p ~ 4e6 + E / (W0 - Wg)
]
reflect = [(Wg - W0) ~ 0] => (terminate_affect!, [Wg], [W0], nothing)
@named bal_eqs = ODESystem(eqs;continuous_events = reflect)
begin
	connections = [
		p1_c.:ψ ~ p1_b.:ψ
		p2_c.:ψ ~ p2_b.:ψ
		p1_b.pressure ~ bal_eqs.p
		p2_b.pressure ~ bal_eqs.p
		bal_eqs.E ~ p1_c.dE + p2_c.dE
		bal_eqs.Wg ~ p1_c.dW + p2_c.dW
	]
	@named connected = ODESystem(connections)
end
composed = compose(
	connected, p1_b, p1_c, p2_b, p2_c, bal_eqs
)
composed_simp = structural_simplify(composed)
u0 = [p1_b.:ψ => 0.0, p1_b.z => 0.0, p2_b.:ψ => 0.0, p2_b.z => 0.0]
pars = [
	p1_b.:κ1=>0.233,
	p1_b.:λ1=>2.26,
	p1_b.:μ1=>0.0,
	p1_b.:κ2=>0.835,
	p1_b.:λ2=>-0.943,
	p1_b.:μ2=>0.0,
	p1_b.Jk=>330e3,
	
	p1_b.zk=>1.53,
	
	p1_c.f=>988e3,
	p1_c.:ρ=>1600,
	p1_c.:α=>1.038e-3,
	p1_c.T_b=>2800,
	p1_c.:ω=>0.05,
	p1_c.:θ=>0.236,

	p2_b.:κ1=>0.233,
	p2_b.:λ1=>2.26,
	p2_b.:μ1=>0.0,
	p2_b.:κ2=>0.835,
	p2_b.:λ2=>-0.943,
	p2_b.:μ2=>0.0,
	p2_b.Jk=>330e3,
	
	p2_b.zk=>1.53,
	
	p2_c.f=>988e3,
	p2_c.:ρ=>1600,
	p2_c.:α=>1.038e-3,
	p2_c.T_b=>2800,
	p2_c.:ω=>0.05,
	p2_c.:θ=>0.236,

	bal_eqs.W0=>0.00125
]
prob = ODEProblem(composed_simp, u0, (0.0, 0.05), pars)
2

I suspect you have not run the code you think you have. As posted, this code does not work since p1_c is not defined. However, accessing W0 works as expected

julia> parameters(bal_eqs)
1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 W0

julia> bal_eqs.W0
bal_eqs₊W0

Working in global scope like this can be pretty confusing, I would restructure this into more functions to make it easier to keep track of what has and has not been executed.

As a side note, converting variable names to symbols is not necessary, p2_b.ω returns the same value as p2_b.:ω.

3

Oh, sorry me, i just forgot to copy on cell from my Pluto notebook. I edit the code

4

OK, I see what is going on. Events are compiled after composition and simplification and only support conditions on state variables of the final system, not observed variables. There is an issue if you want the details.

In your case, adding this to the definition of Wg

begin
    @variables p(t) E(t) Wg(t) [irreducible=true]
    @parameters W0
end

And then adding an initial value to for Wg to u0

u0 = [p1_b.ψ => 0.0, p1_b.z => 0.0, p2_b.ψ => 0.0, p2_b.z => 0.0, bal_eqs.Wg=>0.0]

Allows problem creation to succeed.

5

Thank you a lot!