How can I compute the Jacobian of the system of ODEs obtained by spatially discretising a linear advection equation? A MWE follows:
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
@independent_variables t x
@variables T(..)
Dₜ = Differential(t)
Dₓ = Differential(x)
eq = Dₜ(T(t, x)) ~ -Dₓ(T(t, x))
T₀ = 10.0
bcs = [T(0, x) ~ T₀,
T(t, 0) ~ T₀+exp(-((t-1.0)^2)/0.1)]
domains = [t ∈ Interval(0.0, 1.0),
x ∈ Interval(0.0, 2.0)]
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [T(t, x)])
Δx = 0.1
discretization = MOLFiniteDifference([x => Δx], t, advection_scheme = UpwindScheme())
prob = discretize(pdesys,discretization)
How do I get the Jacobian matrix from the prob variable, which is of type ODEProblem?
An identical question has been asked and answered some two years ago: MethodOfLines.jl: Obtain Jacobian of MOL semi-discretization. The advise given by @ChrisRackauckas there was “During discretize if you pass jac=true then it will generate it, and sys.jac or prob.f.jac are these pieces.” But when I try
prob = discretize(pdesys,discretization,jac=true)
Julia seems to never stop computing. I confess that after some fifteen minutes I had to kill it. The more so that I cannot find such option anywhere in the docs. In fact, I cannot find any docs for the discretize function at all.
Combined with the fact that the advice is now some two years old (maybe it is now deprecated?), that is also why I decided not to “necropost” but rather write a new post.
but Jacobians probably take a lot longer than the right-hand-side functions themselves.
. SymbolicUtils v4 is making Symbolics fully type-stable under the hood, and then we are changing to allow for a good array-based code gen (that’s the end of October goal). Once that is done we will need to update MethodOfLines to hit the array codegen correctly, and then it should have O(1) codegen. Then we really need implement array calculus… otherwise for this exact thing we will need to scalarize. That doesn’t really have a timeline… but anyways that’s how all of this is evolving. So yes the current form is slow because it does it fully on scalars but we have most of the infra to change that now.