Incorporating Auxiliary Calculations as Functions in ModelingToolkit PDE Definitions, Method of Lines

Specific Domains Modelling & Simulations
, , , ,
1

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

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

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.

1 Like