.. |fenics| replace:: FEniCS

.. |petsc| replace:: PETSc

.. |slepc| replace:: SLEPc


.. _solver-choice:

How to choose an appropriate solver?
====================================

|fenics| solver defaults are overzealous, meaning they are biased towards simplicity and robustness (in detriment of performance). A careful dive in the choice of the solver can then have a huge impact in the code performance. This tutorial provides rules of thumb that may be helpful to give a first boost in your solver performance (you may have to dig deeper to get the optimal choice for your specific case).


Parameters are passed to the solver via a Python ``dict`` as:

* a function argument if ``solve`` function is used (``solver_parameters`` argument)
* a class attribute if a solver object (e.g. ``dolfin.cpp.fem.LinearVariationalSolver``) is used (``parameters``)


A simple example for the first case is:


.. code-block:: python

  solver_parameters = {
    'linear_solver': 'cg',
    'preconditioner': 'sor',
  }

  solve(a == L, bcs, solver_parameters=solver_parameters)


For the second case:

.. code-block:: python

  from dolfin.fem.problem import LinearVariationalProblem
  from dolfin.cpp.fem import LinearVariationalSolver

  solver_parameters = {
    'linear_solver': 'cg',
    'preconditioner': 'sor',
  }


  problem = LinearVariationalProblem(a, L, u, bcs)
  solver = LinearVariationalSolver(problem)

  solver.parameters.update(solver_parameters)



Linear problems
---------------

For linear problems (:math:`A x = b`), we lean towards |petsc|. When using this linear algebra backend, the following linear solvers and preconditioners are available, respectively (check out :ref:`intro` for the methods that provide this information):


 ================= =============================================== 
  Name              Method                                         
 ================= =============================================== 
  'bicgstab'        Biconjugate gradient stabilized method         
  'cg'              Conjugate gradient method                      
  'gmres'           Generalized minimal residual method            
  'minres'          Minimal residual method                        
  'petsc'           PETSc built in LU solver                       
  'richardson'      Richardson method                              
  'superlu_dist'    Parallel SuperLU                               
  'tfqmr'           Transpose-free quasi-minimal residual method   
  'umfpack'         UMFPACK                                        
 ================= =============================================== 

|


 ============== ==================================== 
  Name           Method                              
 ============== ==================================== 
  'icc'          Incomplete Cholesky factorization   
  'ilu'          Incomplete LU factorization         
  'petsc_amg'    PETSc algebraic multigrid           
  'sor'          Successive over-relaxation          
 ============== ==================================== 




By default, sparse LU decomposition is used, which is only recommended for systems up to a few thousand unknowns.


To choose the default iterative method, ``linear_solver`` must be set to ``iterative`` or ``krylov``. |fenics| uses `conjugate gradient <https://en.wikipedia.org/wiki/Conjugate_gradient_method>`_ for symmetric linear problems and `generalized minimum residual <https://en.wikipedia.org/wiki/Generalized_minimal_residual_method>`_ for non-symmetric (conjugate gradient is much faster for large problems... that's why we rant so much about :ref:`keeping symmetry <basic-bcs>`). Iterative methods may be sensitive to some parameters (e.g. tolerance).



.. note::

    Check out `this tutorial <https://fenicsproject.org/pub/tutorial/html/._ftut1017.html>`_ for additional information.



Eigenvalue problems
-------------------

For eigenvalue problems, we rely on |slepc|. The possible arguments for ``solver`` are: ``power``, ``subspace``, ``arnoldi``, ``lanczos``, ``krylov-schur``, ``lapack``, ``arpack``, ``jacobi-davidson``, ``generalized-davidson``. You can find `here <https://slepc.upv.es/slepc4py-current/docs/apiref/slepc4py.SLEPc.EPS.Type-class.html>`_ further information about each option.

In :ref:`this tutorial <helmholtz>` a performance comparison between different solvers is done. This may help in your decision process (spoiler: ``generalized-davidson`` is very fast!). 


Other parameters play also a crucial role in performance and convergence ability of the solver. In particular:

* ``problem_type`` should usually be set to ``gen_hermitian`` (the set of available solvers for non-symmetric problems is much smaller - check e.g. `Tzounas, 2020 <https://www.mdpi.com/2076-3417/10/21/7592>`_)

* ``tolerance`` hugely impacts performance


Additionally, ``spectrum`` allows to choose the spectrum region where to look for eigenvalues (e.g. ``smallest magnitude``).