MethodOfLines.jl: Obtain Jacobian of MOL semi-discretization

Specific Domains Modelling & Simulations
1

Consider an exemplary PDE with temporal & spatial derivatives

\partial_t u(t, x) + D_x f\big(u(t, x)\big) = 0

Where D_x is a general differential operator, like the Laplacian and f is some function.

I can use MethodOfLines.jl to obtain a semi-discretization
U(t) = F\big(U(t)\big)
through invoking discretize(sys::PDESystem, disc::D) which gives a AbstractSciMLProblem. Is it possible to call a routine from one of the other SiML packages to obtain the Jacobian
J(U) = \frac{\partial F\big(U(t)\big)}{\partial U} ?
Ideally, this would be possible via algorithmic/automated differentiation, although Finite DIfferences would also be a suitable starting point.

2

During discretize if you pass jac=true then it will generate it, and sys.jac[] or prob.f.jac are these pieces. @xtalax this seems to be missing in the tutorials?

1 Like
3

That is working, thanks!

prob.f.jac returns

(::ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x35fe4ab8, 0xfab08b81, 0x0b2eb902, 0x8f534a9e, 0xca1ca74d)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x00fc76c7, 0xcbe1a643, 0x66f66510, 0x06b3096d, 0xc25da162)}}) (generic function with 2 methods)

Now I would like to evaluate that thing to obtain an actual matrix. I guess one of the arguments arg1, arg2 is the point of evaluation U^* and t is probably time t if the Jacobian is time-dependent, i.e., J\big(t, U(t) \big).
Question is what arg1, arg2 are - I will see if I can figure this out via the ModelingToolkit documentation.

And yeah, in the MethodOfLines.jl repo you do not really find use cases for jac.

4

It’s just the function for DiffEq, so it matches what’s described in the docs.

And because it’s an RGF, you can check its definition via prob.f.jac.body.

1 Like
5

Alright, got it!

One thing I noted is that the evaluation of the jacobian is extremely slow, probably too slow for practical usage.

Although not important for me, prob.f.jac.body returns

ERROR: type #_jac#491 has no field body
Stacktrace:

 [1] getproperty(x::Function, f::Symbol)
   @ Base ./Base.jl:38
6

Its evaluation isn’t slow. It just had bad compilation scaling, which we are working on.