Open ChrisRackauckas opened 4 years ago
Hi, just looked for an open issue where this fits... Is there a way to calculate f
and J
at once ?
I put together my system matrix and the residual vector by a number of calls to ForwardDiff.jacobian
,
so this is coming out naturally. The following works with VoronoiFVM.jl on unstructured grids (I just show the pattern here, will release that code soon):
using DifferentialEquations
using LinearAlgebra
using SparseArrays
mutable struct SolverContext
J
f
uhash
end
function solvercontext(u)
n=length(u)
SolverContext(spzeros(n,n), zeros(n),0x0)
end
function lorenz!(J,du,u)
du[1] = 10.0*(u[2]-u[1])
du[2] = u[1]*(28.0-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)*u[3]
J[1,1]=-10
J[1,2]=10
J[2,1]=28-u[3]
J[2,2]=-1
J[2,3]=-u[1]
J[3,1]=u[2]
J[3,2]=u[1]
J[3,3]=-8/3
end
function lorenz_eval!(p,u)
uhash=hash(u)
if uhash!=p.uhash
lorenz!(p.J,p.f,u)
p.uhash=uhash
end
end
function lorenz_f!(du, u, p,t)
lorenz_eval!(p,u)
du.=p.f
nothing
end
function lorenz_jac!(J, u, p,t)
lorenz_eval!(p,u)
J.=p.J
nothing
end
function run()
u0 = [1.0;0.0;0.0]
ctx=solvercontext(u0)
tspan = (0.0,100.0)
f = ODEFunction(lorenz_f!, jac=lorenz_jac!,jac_prototype=ctx.J)
prob = ODEProblem(f,u0,tspan,ctx)
sol = solve(prob,Rodas5(),reltol=1.0e-10,abstol=1.0e-10)
end
The main idea is that in order to prevent to run the assembly ode twice I just check if it was already run with a given vector by comparing hashes.
Some questions:
Nevertheless I am just lightning struck that everything appears to work: It seems that I can run now my pde evolution problems with DifferentialEquations.jl by just adding a couple of methods to my code which provide the right interface.
And yes, I avoid costly sparsity detection on larger grids here, and my ExtendableSparse.jl matrix class seems to play well, too.
Here is the first running code of this kind (see run_diffeq()) :https://j-fu.github.io/VoronoiFVM.jl/dev/examples/DiffEqCompare/
Is there a way to calculate f and J at once ?
Implicit Euler is probably the only place where that would make sense, and even then I'm a little skeptical. If it takes 20 f
calls to calculate the columns of a matrix, then you at most get <5% speedup, and a Jacobian with 20 columns isn't too uncommon. This interface change would allow you to do these kinds of optimizations though, but this is one that I don't think will make much of a difference because just from the statistics, most of the primal calculations for Jacobians are repeated many times for the columns so the maximal total gain is low.
And yes, I avoid costly sparsity detection on larger grids here, and my ExtendableSparse.jl matrix class seems to play well, too.
Interesting. It might be good to hardcode some matrix coloring solutions for FVM, because indeed it can take a bit to compute. There's an interface for matrix types to overload and provide a fast coloring solution:
https://github.com/SciML/ArrayInterface.jl/blob/master/src/ArrayInterface.jl#L440-L474
How many copies of the Jacobi matrix are around ? IMHO it should be possible to go with just one.
You need two if you're not doing a Krylov method, one if you are. The reason is because factorizations are disconnected from J
updates, so J
and fact(J)
can be at different type points as a way to reduce the total factorizations. This is one of the main ways that the number of factorizations is reduced many fold, so it's a pretty big performance optimization for a 2-fold memory increase. It's common in other stiff ODE solvers as well.
Here is the first running code of this kind (see run_diffeq()) :https://j-fu.github.io/VoronoiFVM.jl/dev/examples/DiffEqCompare/
Thanks for sharing. I think in that case you Rodas5 won't scale well (because it cannot reuse factorizations well). You might want an interface to DAEProblem to use IDA. On the pure Julia side, to really perfect our suite for this kind of problem, RadauIIA5 needs sparse Jacobian support, and we need a pure Julia BDF that is better than the one we had before. We're almost there.
With the new advanced ODE solving tutorial, one thing that has become more clear is that we need to think about how we're doing some of the linear solver pieces.
jac
function that would be approximate, and then use aJacVecOperator
.JacVecOperator
was actually automatic when doing something Newton-Krylov.