Use optimal Newton solver
Many solving operations are currently interfacing with our NewtonSolver class, that builds upon scipy.optimize.newton_krylov
.
This algorithm should be quite fine for sparse matrices (which is the case only for FEM and FD, not for PS discretizations), but I have never checked in detail what other (third-party) solvers we could interface with.
We should checkout different solvers and their performance, e.g., scipy.optimize.root
and the different methods that it implements. In my recent tests, it has been throwing errors, unfortunately.
Also, there is most likely not the optimal solver, but the solving method could differ based on the problem/equation.
Edited by Simon Hartmann