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.
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.
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)
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.
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.
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.
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()