ModelingToolkit model with structural_simplify error: How to initialize variables correctly?

Specific Domains Modelling & Simulations
3

Your problem can be fixed by supplying the defaults to the ODESystem before creating the ODEProblem. Then the only missing variables will be those generated from structural_simplify. I just re-arranged your code a little to supply the defaults like so…

using ModelingToolkit, DifferentialEquations, ForwardDiff, Plots

#Parameters
@parameters L ϵ D μ ρ_g dp
@parameters Rg k1ref Ea1
#variables
@variables z Ca(z) Cb(z) Cc(z) Ctot(z)
@variables y_A(z) y_B(z) y_C(z)
@variables uz(z) T(z) P(z)
@variables R_A(z) R_B(z) R_C(z) Rtotal(z) R1(z)
@variables k1(z)


#Differential Operators
Dz = Differential(z)
#Dzz = Differential(z)^2.0

eqs = [
#Ideal Gas Law
P ~ Ctot * Rg * T
Ca ~ y_A * Ctot
Cb ~ y_B * Ctot 
Cc ~ y_C * Ctot 

#Component mass balances
# 1.0/L * Dz(Ca.*uz) - D*ϵ/(L^2.0) * Dz(Ctot.*Dz(y_A)) ~ R_A
# 1.0/L * Dz(Cb.*uz) - D*ϵ/(L^2.0) * Dz(Ctot.*Dz(y_B)) ~ R_B
# 1.0/L * Dz(Cc.*uz) - D*ϵ/(L^2.0) * Dz(Ctot.*Dz(y_C)) ~ R_C

1.0/L * (Ca.*Dz(uz) + uz.*Dz(Ca)) - D*ϵ/(L^2.0) * (Ctot.*Dz(Dz(y_A))+Dz(y_A).*Dz(Ctot)) ~ R_A
1.0/L * (Cb.*Dz(uz) + uz.*Dz(Cb)) - D*ϵ/(L^2.0) * (Ctot.*Dz(Dz(y_B))+Dz(y_B).*Dz(Ctot)) ~ R_B
1.0/L * (Cc.*Dz(uz) + uz.*Dz(Cc)) - D*ϵ/(L^2.0) * (Ctot.*Dz(Dz(y_C))+Dz(y_C).*Dz(Ctot)) ~ R_C

#Total Mass balance
#1.0/L * Dz(Ctot*uz) ~ Rtotal 
1.0/L * (Ctot*Dz(uz)+ uz.*Dz(Ctot)) ~ Rtotal 

#Pressure Drop
1.0/L * Dz(P) ~  - ρ_g * (1-ϵ)/ϵ^3.0 * uz^2.0/(2*dp)

#Energy balance
#Currently assume isothermal
Dz(T) ~ 0

#Kinetics Equations 
# Arrhenius
k1 ~ k1ref * exp(-Ea1./Rg./T)

#rates
R1 ~ (k1.*(Ca^0.5))/(1.0 + k1.*Ca^0.5)

#Molar reaction rates
R_A ~ -1.0*R1 
R_B ~ -0.5*R1
R_C ~ +1.0*R1

#Total reaction rate
Rtotal ~ R_A + R_B + R_C 
]

bcs =[
	  Ca => 1.74,
	  Cb => 34.0,
	  Cc => 0.0,
      Ctot => 34.0+1.74,
	  T => 517.0,
      P => (34.0+1.74)*Rg*517.0,
      y_A => 1.74/(34.0+1.74),
      y_B => 34.0/(34.0+1.74),
      y_C => 0.0,
	  k1 => k1ref.*exp(-Ea1/Rg/517.0),
	  R1 => 1E-10,
      uz => 0.47,
      R_A => 0.0,
      R_B => 0.0,
      R_C => 0.0,
      Rtotal => 0.0,
      Dz(y_A) => 0.0,
      Dz(y_B) => 0.0,
      Dz(y_C) => 0.0,
	  ]

#Parameter map
pars = [
Rg => 8.314, 
ϵ  => 0.4,
Ea1 => 79.496*1000,
k1ref => 125.0E7*0.2,
L  => 0.7,
dp => 0.0046,
ρ_g  => 1.018,
D  => 4.9E-6,
μ => 1E-8
] 

@named pfr = ODESystem(eqs, z; defaults=[bcs; pars])

Now you can properly assembly your problem

Lspan = (0.0,1.0)
simpsys = structural_simplify(pfr)

missing_bcs = ModelingToolkit.missing_variable_defaults(simpsys)

prob = ODEProblem(simpsys,missing_bcs,Lspan,pars,jac=true)

The Rodas solver didn’t work for me. I was only able to get it to solve with the non-adaptive ImplicitEuler, or you can use the SimpleEuler* package to get a higher order solve…

sol = solve(prob, ImplicitEuler(nlsolve=NLNewton(check_div=false)); dt=1e-3, adaptive=false, isoutofdomain=(u,p,t)->any(x->x<0,u))

using SimpleEuler
sol = solve(prob, BackwardEuler(order=2); isoutofdomain=(u,p,t)->any(x->x<0,u))
plot(sol; idxs=[y_A, y_C])

image