$$\newcommand{\dt}{\Delta t} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\x}{\boldsymbol{x}} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\uI}{u_{_0}} \newcommand{\ub}{u_{_\mathrm{D}}} \newcommand{\GD}{\Gamma_{_\mathrm{D}}} \newcommand{\GN}{\Gamma_{_\mathrm{N}}} \newcommand{\GR}{\Gamma_{_\mathrm{R}}} \newcommand{\inner}{\langle #1, #2 \rangle}$$

# Working with linear solvers

Sparse LU decomposition (Gaussian elimination) is used by default to solve linear systems of equations in FEniCS programs. This is a very robust and simple method. It is the recommended method for systems with up to a few thousand unknowns and may hence be the method of choice for many 2D and smaller 3D problems. However, sparse LU decomposition becomes slow and one quickly runs out of memory for larger problems. For large problems, we instead need to use iterative methods which are faster and require much less memory. We will now look at how to take advantage of state-of-the-art iterative solution methods in FEniCS.

## Choosing a linear solver and preconditioner

Preconditioned Krylov solvers is a type of popular iterative methods that are easily accessible in FEniCS programs. The Poisson equation results in a symmetric, positive definite system matrix, for which the optimal Krylov solver is the Conjugate Gradient (CG) method. For non-symmetric problems, a Krylov solver for non-symmetric systems, such as GMRES, is a better choice. Incomplete LU factorization (ILU) is a popular and robust all-round preconditioner, so let us try the GMRES-ILU pair:

solve(a == L, u, bc,
solver_parameters={'linear_solver': 'gmres',
'preconditioner': 'ilu'})
# Alternative syntax
solve(a == L, u, bc,
solver_parameters=dict(linear_solver='gmres',
preconditioner='ilu'))


the section List of linear solver methods and preconditioners lists the most popular choices of Krylov solvers and preconditioners available in FEniCS.

## Choosing a linear algebra backend

The actual GMRES and ILU implementations that are brought into action depend on the choice of linear algebra package. FEniCS interfaces several linear algebra packages, called linear algebra backends in FEniCS terminology. PETSc is the default choice if FEniCS is compiled with PETSc. If PETSc is not available, then FEniCS falls back to using the Eigen backend. The linear algebra backend in FEniCS can be set using the following command:

parameters.linear_algebra_backend = backendname


where backendname is a string. To see which linear algebra backends are available, you can call the FEniCS function list_linear_algebra_backends. Similarly, one may check which linear algebra backend is currently being used by the following command:

print(parameters.linear_algebra_backend)


## Setting solver parameters

We will normally want to control the tolerance in the stopping criterion and the maximum number of iterations when running an iterative method. Such parameters can be controlled at both a global and a local level. We will start by looking at how to set global parameters. For more advanced programs, one may want to use a number of different linear solvers and set different tolerances and other parameters. Then it becomes important to control the parameters at a local level. We will return to this issue in the section Linear variational problem and solver objects.

Changing a parameter in the global FEniCS parameter database affects all linear solvers (created after the parameter has been set). The global FEniCS parameter database is simply called parameters and it behaves as a nested dictionary. Write

info(parameters, verbose=True)


to list all parameters and their default values in the database. The nesting of parameter sets is indicated through indentation in the output from info. According to this output, the relevant parameter set is named 'krylov_solver', and the parameters are set like this:

prm = parameters.krylov_solver  # short form
prm.absolute_tolerance = 1E-10
prm.relative_tolerance = 1E-6
prm.maximum_iterations = 1000


Stopping criteria for Krylov solvers usually involve some norm of the residual, which must be smaller than the absolute tolerance parameter or smaller than the relative tolerance parameter times the initial residual.

We remark that default values for the global parameter database can be defined in an XML file. To generate such a file from the current set of parameters in a program, run

File('parameters.xml') << parameters


If a dolfin_parameters.xml file is found in the directory where a FEniCS program is run, this file is read and used to initialize the parameters object. Otherwise, the file .config/fenics/dolfin_parameters.xml in the user's home directory is read, if it exists. Another alternative is to load the XML file (with any name) manually in the program:

File('parameters.xml') >> parameters


The XML file can also be in gzip'ed form with the extension .xml.gz.

## An extended solver function

We may extend the previous solver function from ft12_poisson_solver.py in the section A more general solver function such that it also offers the GMRES+ILU preconditioned Krylov solver:

This new solver function, found in the file ft10_poisson_extended.py, replaces the one in ft12_poisson_solver.py. It has all the functionality of the previous solver function, but can also solve the linear system with iterative methods.

## A remark regarding unit tests

Regarding verification of the new solver function in terms of unit tests, it turns out that unit testing for a problem where the approximation error vanishes gets more complicated when we use iterative methods. The problem is to keep the error due to iterative solution smaller than the tolerance used in the verification tests. First of all, this means that the tolerances used in the Krylov solvers must be smaller than the tolerance used in the assert test, but this is no guarantee to keep the linear solver error this small. For linear elements and small meshes, a tolerance of $$10^{-11}$$ works well in the case of Krylov solvers too (using a tolerance $$10^{-12}$$ in those solvers). The interested reader is referred to the demo_solvers function in ft10_poisson_extended.py for details: this function tests the numerical solution for direct and iterative linear solvers, for different meshes, and different degrees of the polynomials in the finite element basis functions.

## List of linear solver methods and preconditioners

Which linear solvers and preconditioners that are available in FEniCS depends on how FEniCS has been configured and which linear algebra backend is currently active. The following table shows an example of which linear solvers that can be available through FEniCS when the PETSc backend is active:

 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

The set of available preconditioners also depends on configuration and linear algebra backend. The following table shows an example of which preconditioners may be available:

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

An up-to-date list of the available solvers and preconditioners for your FEniCS installation can be produced by

list_linear_solver_methods()
list_krylov_solver_preconditioners()