Hello everyone,

I’ve been working on a PDE model using the `ModelingToolkit`

, `MethodOfLines`

, and `DifferentialEquations`

packages, and I’ve run into some challenges that I’m hoping to get guidance on. For confidentiality reasons, I can’t share the exact nature of the equations I’m working with, but I’ll do my best to describe the problem.

I have a series of PDEs which depend on multiple auxiliary calculations, many of which are inter-related. The direct inclusion of these calculations in my PDE definitions leads to a long and unwieldy code, and I believe it may be affecting performance. I could not figure out how to incorporate these auxiliary calculations. I’ve been advised to incorporate these auxiliary calculations as separate functions to streamline the code, but I’m struggling with how to best do this in the `ModelingToolkit`

framework. Or is there another way to do this?

I’m unsure of how to incorporate these functions into the PDE definitions, especially given that some functions rely on the outputs of others.

- Can someone guide me on best practices for incorporating auxiliary calculations as functions in
`ModelingToolkit`

PDE definitions and usage with the`MethodOfLines.jl`

? - Any recommendations for ensuring that the performance is optimized when using these auxiliary functions?
- When I use all of these auxiliary equations, it leads to instability issues even when I define tolerances
`(abstol=1.0E-8,reltol=1.0E-8)`

with`Rosenbrock23()`

. Is there a way to work around this without reducing the complexity of the model?

Here’s a simplified version of my code to give you a rough idea (I am illustrating only one equation and related auxiliary equations, boundary initial conditions):

```
using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets
@parameters t z z_sub
@variables I(..) Id(..) λ_0(..) Z(..) Vm(..) Is(..) T(..) T_sub(..)
Dt = Differential(t)
Dz = Differential(z)
Dzz = Differential(z)^2
Dz_sub = Differential(z_sub)
Dzz_sub = Differential(z_sub)^2
t_min = z_min = 0.0
z_max = p.film_thickness
t_max = 30.0
z_sub_min = p.film_thickness
z_sub_max = p.film_thickness + p.thickness_substrate
# Auxiliary equations
M = Vm(t,z)/2
vT = M*p.MW_M / p.rho_M + (M0 - M)*p.MW_M / p.rho_P
phi_M = M*p.MW_M / (p.rho_M*vT)
phi_P = (M0 - M)*p.MW_M / (p.rho_P*vT)
f = p.f0 + phi_M * p.α_M * (T(t,z) - p.Tg_M) + phi_P * p.α_P * (T(t,z) - p.Tg_P)
I_light = current_intensity*exp(-ϵ*Id(t,z)*z-ϵ_PIG*c_p*z)
kd = (ϵ*Φ*I_light)/E_photon
D_i = p.Di0*exp(-p.A_Di /f)
kt=...
D_z=...
...
eq = [
Dt(I(t,z)) ~ Dz(D_i*Dz(I(t,z))) - kd*I(t,z),
.... # Other equations
]
bcs = [
I(0,z) ~ 105.2191,
Dz(I(t,0)) ~ 0.0,
Dz(I(t,z_max)) ~ 0.0,
... # Other bc
]
@named pdesys = PDESystem(eq,bcs,domains,[z,z_sub,t],[I(t,z),Id(t,z),λ_0(t,z),Vm(t,z),Z(t,z),Is(t,z),T(t,z),T_sub(t,z_sub)])
discretization = MOLFiniteDifference([z=>50,z_sub=>55],t,approx_order=2)
pred = solve(prob,Rosenbrock23();saveat=0.01,abstol=1.0E-10,reltol=1.0E-10)
```