MethodOfLines.jl: Obtain Jacobian of MOL semi-discretization

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.

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

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.

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

https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/#SciMLBase.ODEFunction

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

1 Like

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

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