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

Hello everybody,
I am currently developing a ModelingToolkit.jl model to simulate axial flow in a tubular pipe with a chemical reaction, including axial diffusion. This system is currently an initial value problem, and I want to simulate the evolution of several variables along the length of pipe. The general molar balance for species j in the fluid mixture is written as (apologize for the screenshots as I do not know how to write equations in this forum):

image

Here, the variables that are in bold fond are functions of axial position in the pipe and the normal typeface denotes parameters. As you can see on the LHS, there are variables inside derivatives. I could not get Julia to automatically handle two variables inside a derivative term, so I used a product rule to expand these terms to the equation below:

image

The RHS term denotes a chemical reaction term A+B → C, which is an algebraic function:

image

The rest of the equations are presented in the simplified model shown at the end of this post. In this current model, I am running into issues with structural_simplify. Once I run the command and try to create a new ODEProblem with the simplified system, I get the following error message:

ERROR: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[Ctotˍz(z), y_Aˍzz(z), y_Bˍzz(z), y_Cˍzz(z), uzˍz(z)] are missing from the variable map.

As far as I can tell, this is being caused due to the new variables structural_simplify generates during simplification. I looked at this section of the ModelingToolkit.jl FAQ, which suggests running the command

ModelingToolkit.missing_variable_defaults()

for my simplified system, unfortunately, this command makes all initial guesses for my system 0.0, even for the initial conditions I set for the original model. Not sure if I am using this incorrectly for the simplified model? The FAQ shows you can hard code a vector of initial guesses instead of zeroes, but depending on a given model the order of these variables and which variables are created may differ. Is there a way to just add the default IC values for the variables ModelingToolkit generates? I have included initial condition values for Dz(y_A), Dz(y_B), and Dz(y_C), as those are the variables that have a double derivative in the mole balance equation, yet it seems this initial condition is insufficient. Any feedback on how to properly specify the initial conditions for a model that uses structural_simplify would be appreciated!

Finally, as a related note, I am wondering about the need for initial conditions for algebraic variables. Is there a way to make ModelingToolkit calculate the initial values of algebraic variables based on the initial condition for differential variables? In this scenario for example, If I know the initial condition for Ca, T, and P, in principle the values for R_j, Ctot, and R_total are known. Is there a best practice to automatically calculate the initial conditions of these variables? Right now, I initialize all algebraic variables with 0.0 or with hard-coded equations, but perhaps there may be instabilities if these initial conditions are too far-off the values that the actual algebraic equations calculate.

Thank you for any feedback or suggestions on how to proceed! The code is shown below:

using ModelingToolkit, DifferentialEquations, NonlinearSolve, 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

@named pfr = ODESystem([
#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
] 

Lspan = (0.0,1.0)
simpsys = structural_simplify(pfr)
prob = ODEProblem(simpsys,bcs,Lspan,pars,jac=true)
# prob = ODEProblem(simpsys,ModelingToolkit.missing_variable_defaults(simpsys),Lspan,pars,jac=true) #ODEProblem with default values for ICs
sol = solve(prob,Rodas5P(),reltol=1E-10,abstol=1E-10,isoutofdomain=(u,p,t)->any(x->x<0,u))

Yes, those are dummy derivative terms. You need initial conditions for D(Ctot) and such in order for the DAE to be well-specified.

This is something I’m addressing in Fix up initializesystem for hierarchical models by ChrisRackauckas · Pull Request #2403 · SciML/ModelingToolkit.jl · GitHub with the ModelingToolkit v9 release that is schedule to be finalizing in the next week or so.

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

1 Like

Thanks for working through the solution! The part I am not following fully in your solution is where simpsys uses the defaults provided. when you set defaults=[bcs; pars], and then use the missing_bcs, does it only set the new default values to any variable that was not in bcs?

Interesting that the ImplicitEuler was the only way you made it work, I will incorporate your changes and attempt that solver as well.

Update: I was able to reproduce the successful solving of the model with ImplicitEuler. I wonder why Rodas solvers don’t work with the adaptive time-stepping.

Thanks for the replies! I got the script to complete the ODEProblem by including the differential initial conditions for all variables that contain a derivative, as well as double derivative initial; conidtions for y_A, y_B, and y_C, which were the terms with double derivatives. Not sure I understand why, as I believed you only needed differential initial conditions for higher-order derivatives, which would only be y_A, y_B, and y_C, and would not need double derivative initial conditions. Would this be suggestive of another issue on how I set up the ModelingToolkit problem?

The initial conditions are now written as:


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,
      Dz(Ctot) => 0.0,
      Dz(Ca) => 0.0,
      Dz(Cb) => 0.0,
      Dz(Cc) => 0.0,
      Dz(uz) => 0.0,
      Dz(P) => 0.0,
      Dz(T) => 0.0,
      Dz(Dz(y_A)) => 0.0,
      Dz(Dz(y_B)) => 0.0,
      Dz(Dz(y_C)) => 0.0,
	  ]

Obviously, its more elegant to use Brad_Carman’s solution until I can test out the ModelingToolkit update you mention for initialization of algebraic equations, which sounds great!

One additional update: The Rodas5P solver works if I set the solver to not use adaptive time-stepping. setting dt=1E-4 or lower leads to a domain error issue:

ERROR: DomainError with -14.058907925955667:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

Is this issue and the adaptive time-stepping issue diagnostic of other model problems? Why may the adaptive algorithms run into problems near the initial condition at L=0? My thought here is that the time step dt = 1E-3 ‘skips’ the region where the integration runs into some sort of problem, allowing the system to continue. Not sure how to troubleshoot what is going on there.