I’m bringing up a challenging but pragmatic question. How to choose the best solver for a given problem? At least, considering the FEniCS framework or its default linear algebra backend, PETSc.
As a beginner, I’m asking this because I often spend much time finding the optimal solver for my problems (given the mesh, the set of PDEs, the parameters, the discretization, and so on) through trial and error. What I learned from basic theory (like this book) is the general rule of thumb:
- direct solvers (e.g.
mumps) → problems on small/medium 2D domains
- iterative solvers (e.g.
gmres) → problems on medium/large 2D and 3D domains
So, given the dimension of the problem, I simply try all the solvers available of the appropriate type (direct/iterative), and I select the fastest among those which converge.
This approach also seems supported by Matthew Knepley (PETSc developer), who in this talk clearly states that you never know which is the best solver for your problem.
Still, I was wondering if this is the most efficient way. So, rephrasing my main questions in more specific (and bearable) subquestions:
- Do you think trial and error is the best way to find the optimal solver for a problem?
- Are there valuable rules of thumb that I can use to select the optimal solver a priori or, at least, select a subset in advance?
- Are there any references you’d recommend to me about this topic?
Thank you very much.