# Incorporating Auxiliary Calculations as Functions in ModelingToolkit PDE Definitions, Method of Lines
**URL:** https://discourse.julialang.org/t/incorporating-auxiliary-calculations-as-functions-in-modelingtoolkit-pde-definitions-method-of-lines/104428
**Category:** Modelling & Simulations
**Tags:** question, pde, modelingtoolkit, differentialequation, methodoflines
**Created:** [September 30, 2023, 12:07pm UTC](https://discourse.julialang.org/t/incorporating-auxiliary-calculations-as-functions-in-modelingtoolkit-pde-definitions-method-of-lines/104428 "2023-09-30T12:07:22Z")
**Posts on this page:** 2
**Page:** 1
### Author:  [@Can\_Bozdogan](https://discourse.julialang.org/u/Can_Bozdogan)
#### Post date: [September 30, 2023, 12:07pm UTC](https://discourse.julialang.org/t/incorporating-auxiliary-calculations-as-functions-in-modelingtoolkit-pde-definitions-method-of-lines/104428/1 "2023-09-30T12:07:22Z")
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.
1. Can someone guide me on best practices for incorporating auxiliary calculations as functions in `ModelingToolkit` PDE definitions and usage with the `MethodOfLines.jl`?
2. Any recommendations for ensuring that the performance is optimized when using these auxiliary functions?
3. 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):
```julia
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)
```
---
### Author:  [@xtalax](https://discourse.julialang.org/u/xtalax)
#### Post date: [October 3, 2023, 1:09pm UTC](https://discourse.julialang.org/t/incorporating-auxiliary-calculations-as-functions-in-modelingtoolkit-pde-definitions-method-of-lines/104428/2 "2023-10-03T13:09:21Z")
You can write a function that takes mtk variables as input and performs arbitrary calculations with them, also leading to optimised code once generated. The only caveat being that it must not contain branches, and if it needs to then you have to use `ifelse` instead. If you can’t do this, you can also register the function with `@register_symbolic f(t, x, u, v...)` which will stop symbolic tracing, but then you need to take care over the efficiency of your code.