Unexpectedly Large Compile Time Using MethodOfLines.jl

Specific Domains Numerics
1

So I am trying to use the MethodOfLines package to solve some equations which are somewhat similar to 2D diffusion equations. As a starting point, I have modified the heat equation example from the documentation to work in 2D, here is my code:

module Diffusion2D

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
# Parameters, variables, and derivatives
@parameters t x y
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
Dy = Differential(y)
Dyy = Differential(y)^2

# 1D PDE and boundary conditions
eq  = Dt(u(t, x, y)) ~ 0.1*(Dxx(u(t, x, y)) + Dyy(u(t, x, y)))
bcs = [u(0, x, y) ~ sin(2*pi*x) * sin(2*pi*y),
        u(t, 0, y) ~ 0.0,
        u(t, 1, y) ~ 0.0,
        u(t, x, 0) ~ 0.0,
        u(t, x, 1) ~ 0.0]
 
# Space and time domains
domains = [t ∈ Interval(0.0, 0.5),
        x ∈ Interval(0.0, 1.0),
        y ∈ Interval(0.0, 1.0)]

# PDE system
@named pdesys = PDESystem(eq, bcs, domains,[t, x, y],[u(t, x, y)])

# Method of lines discretization
dx = 0.01
dy = 0.01
order = 2
discretization = MOLFiniteDifference([x => dx, y => dy],t, approx_order=order)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization)

# Solve ODE problem
@time sol = solve(prob, Tsit5(), saveat=0.01)

end

This code seems to work correctly and it gives the right solution, and once it’s been compiled, the solve step runs in ~140 seconds typically, which is within reason for the dx/dy I am using. The main issue is that the first time I run the code, the solve step typically takes ~6000 seconds to compile and run (it was even longer than this before I put it all in a module). This strikes me as an unusually large compile time, have I made a mistake which is causing this issue? Or otherwise is there ways I can optimize the compile time? If it’s relevant, I am running this code inside a Jupyter notebook. Thanks!

2

The current codegen of MethodOfLines uses a scalarization which leads to long compile times. This is something that is being addressed by JuliaSimCompiler:

We will demonstrate this in the near future.

2 Likes
3

Ah I see, thanks for the response! Is there currently a way to convert PDESystems to IRSystems in the same way that is done with the ODESystem used in the example on the JuliaSimCompiler page (or at least achieve the same effect)? I tried modifying my above code to add:

...

@named pdesys = PDESystem(eq, bcs, domains,[t, x, y],[u(t, x, y)])

using JuliaSimCompiler
pdesys_ir = IRSystem(pdesys)

...

but it seems that IRSystem() does not currently have a method for PDESystem.

1 Like
4

There is an extension being created so that loading JuliaSimCompiler will make it automatically use IRSystem. These are updates that will be demonstrated in the near future.