[ANN] MethodOfLines v1: Symbolic Arrays => Massive Compile Speedups

Package Announcements
, ,
1

MethodOfLines.jl v1.0 changes all discretizations to using symbolic arrays. If they lower to code using DAEProblem the array structure is kept intact, meaning no more scalarization and massive compile time speedups.

Timings

Earlier benchmarks comparing the pointwise symbolic representation with the new array form measured the following discretization times:

Problem Unknowns Pointwise form Array form Speedup
1D heat equation 400 47.9 s 0.23 s 201×
2D heat equation 900 7.52 s 0.14 s 52×
3D heat equation 1,331 7.78 s 0.13 s 60×

The same effect applies to nonlinear operators. For

Dt(u(t, x)) ~ Dx(u(t, x) * Dx(u(t, x)))

the results were:

Grid points Pointwise expression nodes Array expression nodes Pointwise time Array time Speedup
21 1,653 376 0.04 s 0.02 s 2×
101 8,613 376 1.20 s 0.05 s 24×
501 43,413 376 49.4 s 0.44 s 112×

The numerical solutions agreed to within 9e-16.

I also measured the complete problem-construction path on the current v1 code. These are single warmed runs on Julia 1.12.6 using an AMD EPYC 7502 system:

Grid points New DAE path Compiled ODE path Speedup
256 0.59 s, 61 MiB 1.56 s, 165 MiB 2.7×
1,024 1.77 s, 315 MiB 12.58 s, 1.31 GiB 7.1×
4,096 17.48 s, 5.90 GiB 175.22 s, 20.63 GiB 10.0×

Caveat: Using explicit Runge–Kutta methods

The array compilation advantage currently applies to the DAE path. Explicit Runge–Kutta algorithms such as Tsit5() and SSPRK54() require an ODEProblem.

To use an explicit method, request the symbolic system directly and compile it into an ODE:

sys, tspan = symbolic_discretize(pdesys, disc)
prob = ODEProblem(mtkcompile(sys), nothing, tspan)

sol = solve(prob, Tsit5())

This path scalarizes the system during mtkcompile, so it does not have the same constant-size symbolic representation or compilation scaling as the default DAE path. We plan to address this limitation in the near future.

2

A kind of naive question… I am working on some problems where I make components consisting of PDEs (e.g., gas pipes). In my set-up, I need to be able to have ports for the pipes, and be able to connect several pipes to some junction. For some problems, I have gas mixtures and need to invoke thermodynamic models.

Is it feasible to use MOL for such pipes? I have currently done some discretization on my own (with AI tools…), with both method of characteristic (MOC) and FVM/Kurganov-Petrova 2nd order method. I found it really tricky to get this to work, partially because I am not an expert on numeric methods for PDEs.

3

I haven’t tried to use MethodOfLines.jl in a component, but that and UDEs in MOL are two next direction

4

Two experiences I had with implementing PDE discretization and one general experience…

  1. One tricky part I found was to choose proper boundary conditions that made the port/connector work.
  2. I also struggled quite a bit with making the problem well scaled. I like to write PDEs in the “original”/“balance equation” variables, so I prefer to use mass per length as differential variable (vector) in a 1D formulation of the mass balance, and momentum per length. Then add AEs. Reason: I find it easier to debug the model. Google AI Studio insisted I should use pressure and velocity (?) because it thought that made the problem better posed. (I have started to experiment with Claude now).
  3. A recurrent problem I have had with MTK, and even more when I started to work with components, is that of properly initializing the problem … to avoid over specified or under specified initialization. I think I have found a solution to this, and I also think this solution (or a similar one) would be very useful to MTK users in general, see example below.
u0_map, guess_map = build_mtk_maps(csys, sys)
update_mtk_map(u0_map)
update_mtk_map(guess_map)
prob = ODEProblem(csys, u0_map, tspan; guesses=guess_map)

assumes the compiled system csys is index 1 DAE. build_mtk_maps pulls out unknowns from csys and uses equations to check which unknown is differential. The differential variables u0_map are then set equal to the specified initial condition or guess value (initial condition is prioritized), and all algebraic variables are set to nothing. Similar for the guess_map wrt. initial conditions/guess values (nothing is of course not introduced here). Second argument sys is used to detect extra variables introduced in the index reduction process, and these are set to 0.0 for lack of something better.

If some of the unknowns lack initial conditions or initial guesses, they are set to missing.

Function update_mtk_map with one argument lists variables that have value missing– the user needs to specify these; this is done by adding a map as a the second argument of update_mtk_map.

Anyways, I find these functions highly useful, and since I developed them (with Google AI Studio, and improved with Claude), there is no more problem with overdetermined or underdetermined initial conditions. I also have a twin build_mtk_maps which uses a solution as the second argument instead of sys – this is used to re-start the system from the final value of a previous solution.

The only problem that remains is to select good initial guesss for algebraic variables…

5

If I just create an MTK component, not using MethodOfLines at all, and write my equation in array form, will they still be automatically scalarized (this is what was happening if I understand correctly) or will they be handled as array equation?

6

If using DAEProblem and no index reduction, then it will not scalarize