Is it possible to solve a nonlinear problem matrix-free?

Question

I have a nonlinear problem where I can easily compute F(u), but its Jacobian is dense.

Since computing each Jacobian-vector product is cheap, I'd like to wire this up (with NonlinearProblem/NewtonSolver?) to solve the nonlinear problem with a matrix-free Krylov iteration.

(There used to be some partial support for matrix-free in DOLFIN 1.0, but I think it was removed.)

Any suggestions or comments gratefully received.

If it isn't supported (which I expect), can we talk about its potential design? (I'm sure Garth, Anders etc will have opinions.) Maybe instead of

problem = MyNonlinearProblemSubclass()
solver = NewtonSovler("cg")
solver.solve(problem, u.vector()) # calls problem.J to compute sparse Jacobian

we could have

problem = MyNonlinearProblemSubclass()
solver = NewtonSolver("cg", matrix_free=True)
solver.solve(problem, u.vector()) # calls problem.Jv to compute Jacobian action

?

Patrick

Answer 1

There is no matrix-free demo, but there are some unit tests in test/unit/la/python/LinearOperator.py that illustrate a matrix-free Krylov linear solver. I don't know how well it will hook up to the NewtonSolver.

I have some thoughts on how one might want to pass a matrix-free solver into a NewtonSolver, but it would be best to discuss this on the mailing list.