# ModelingToolkit KeyError

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)
``````

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.:ω`.

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

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.

Thank you a lot!