$$\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}$$

# High-level and low-level solver interfaces

The FEniCS interface allows different ways to access the core functionality, ranging from very high-level to low-level access. So far, we have mostly used the high-level call solve(a == L, u, bc) to solve a variational problem a == L with a certain boundary condition bc. However, sometimes you may need more fine-grained control of the solution process. In particular, the call to solve will create certain objects that are thrown away after the solution has been computed, and it may be practical or efficient to reuse those objects.

## Linear variational problem and solver objects

In this section, we will look at an alternative interface to solving linear variational problems in FEniCS, which may be preferable in many situations compared to the high-level solve function interface. This interface uses the two classes LinearVariationalProblem and LinearVariationalSolver. Using this interface, the equivalent of solve(a == L, u, bc) looks as follows:

u = Function(V)
problem = LinearVariationalProblem(a, L, u, bc)
solver = LinearVariationalSolver(problem)
solver.solve()


Many FEniCS objects have an attribute parameters, similar to the global parameters database, but local to the object. Here, solver.parameters play that role. Setting the CG method with ILU preconditioning as the solution method and specifying solver-specific parameters can be done like this:

solver.parameters.linear_solver = 'gmres'
solver.parameters.preconditioner = 'ilu'
prm = solver.parameters.krylov_solver  # short form
prm.absolute_tolerance = 1E-7
prm.relative_tolerance = 1E-4
prm.maximum_iterations = 1000


Settings in the global parameters database are propagated to parameter sets in individual objects, with the possibility of being overwritten as above. Note that global parameter values can only affect local parameter values if set before the time of creation of the local object. Thus, changing the value of the tolerance in the global parameter database will not affect the parameters for already created solvers.

## Explicit assembly and solve

As we saw already in the section The Navier--Stokes equations, linear variational problems can be assembled explicitly in FEniCS into matrices and vectors using the assemble function. This allows even more fine-grained control of the solution process compared to using the high-level solve function or using the classes LinearVariationalProblem and LinearVariationalSolver. We will now look more closely into how to use the assemble function and how to combine this with low-level calls for solving the assembled linear systems.

Given a variational problem $$a(u,v)=L(v)$$, the discrete solution $$u$$ is computed by inserting $$u=\sum_{j=1}^N U_j \phi_j$$ into $$a(u,v)$$ and demanding $$a(u,v)=L(v)$$ to be fulfilled for $$N$$ test functions $$\hat\phi_1,\ldots,\hat\phi_N$$. This implies $$\begin{equation*} \sum_{j=1}^N a(\phi_j,\hat\phi_i) U_j = L(\hat\phi_i),\quad i=1,\ldots,N, \end{equation*}$$ which is nothing but a linear system, $$\begin{equation*} AU = b, \end{equation*}$$ where the entries of $$A$$ and $$b$$ are given by \begin{align*} A_{ij} &= a(\phi_j, \hat{\phi}_i), \\ b_i &= L(\hat\phi_i)\tp \end{align*}

The examples so far have specified the left- and right-hand sides of the variational formulation and then asked FEniCS to assemble the linear system and solve it. An alternative is to explicitly call functions for assembling the coefficient matrix $$A$$ and the right-hand side vector $$b$$, and then solve the linear system $$AU=b$$ for the vector $$U$$. Instead of solve(a == L, U, b) we now write

A = assemble(a)
b = assemble(L)
bc.apply(A, b)
u = Function(V)
U = u.vector()
solve(A, U, b)


The variables a and L are the same as before; that is, a refers to the bilinear form involving a TrialFunction object u and a TestFunction object v, and L involves the same TestFunction object v. From a and L, the assemble function can compute $$A$$ and $$b$$.

Creating the linear system explicitly in a program can have some advantages in more advanced problem settings. For example, $$A$$ may be constant throughout a time-dependent simulation, so we can avoid recalculating $$A$$ at every time level and save a significant amount of simulation time.

The matrix $$A$$ and vector $$b$$ are first assembled without incorporating essential (Dirichlet) boundary conditions. Thereafter, the call bc.apply(A, b) performs the necessary modifications of the linear system such that u is guaranteed to equal the prescribed boundary values. When we have multiple Dirichlet conditions stored in a list bcs, we must apply each condition in bcs to the system:

for bc in bcs:
bc.apply(A, b)

# Alternative syntax using list comprehension
[bc.apply(A, b) for bc in bcs]


Alternatively, we can use the function assemble_system, which takes the boundary conditions into account when assembling the linear system. This method preserves the symmetry of the linear system for a symmetric bilinear form. Even if the matrix A that comes out of the call to assemble is symmetric (for a symmetric bilinear form a), the call to bc.apply will break the symmetry. Preserving the symmetry of a variational problem is important when using particular linear solvers designed for symmetric systems, such as the conjugate gradient method.

Once the linear system has been assembled, we need to compute the solution $$U=A^{-1}b$$ and store the solution $$U$$ in the vector U = u.vector(). In the same way as linear variational problems can be programmed using different interfaces in FEniCS—the high-level solve function, the class LinearVariationalSolver, and the low-level assemble function—linear systems can also be programmed using different interfaces in FEniCS. The high-level interface to solving a linear system in FEniCS is also named solve:

solve(A, U, b)


By default, solve(A, U, b) uses sparse LU decomposition to compute the solution. Specification of an iterative solver and preconditioner can be made through two optional arguments:

solve(A, U, b, 'cg', 'ilu')


Appropriate names of solvers and preconditioners are found in the section List of linear solver methods and preconditioners.

This high-level interface is useful for many applications, but sometimes more fine-grained control is needed. One can then create one or more KrylovSolver objects that are then used to solve linear systems. Each different solver object can have its own set of parameters and selection of iterative method and preconditioner. Here is an example:

solver = KrylovSolver('cg', 'ilu')
prm = solver.parameters
prm.absolute_tolerance = 1E-7
prm.relative_tolerance = 1E-4
prm.maximum_iterations = 1000
u = Function(V)
U = u.vector()
solver.solve(A, U, b)


The function solver_linalg in the program file ft10_poisson_extended.py implements such a solver.

The choice of start vector for the iterations in a linear solver is often important. By default, the values of u and thus the vector U = u.vector() will be initialized to zero. If we instead wanted to initialize U with random numbers in the interval $$[-100,100]$$ this can be done as follows:

n = u.vector().array().size
U = u.vector()
U[:] = numpy.random.uniform(-100, 100, n)
solver.parameters.nonzero_initial_guess = True
solver.solve(A, U, b)


Note that we must both turn off the default behavior of setting the start vector ("initial guess") to zero, and also set the values of the vector U to nonzero values. This is useful if we happen to know a good initial guess for the solution.

Using a nonzero initial guess can be particularly important for time-dependent problems or when solving a linear system as part of a nonlinear iteration, since then the previous solution vector U will often be a good initial guess for the solution in the next time step or iteration. In this case, the values in the vector U will naturally be initialized with the previous solution vector (if we just used it to solve a linear system), so the only extra step necessary is to set the parameter nonzero_initial_guess to True.

## Examining matrix and vector values

When calling A = assemble(a) and b = assemble(L), the object A will be of type Matrix, while b and u.vector() are of type Vector. To examine the values, we may convert the matrix and vector data to numpy arrays by calling the array method as shown before. For example, if you wonder how essential boundary conditions are incorporated into linear systems, you can print out A and b before and after the bc.apply(A, b) call:

A = assemble(a)
b = assemble(L)
if mesh.num_cells() < 16:  # print for small meshes only
print(A.array())
print(b.array())
bc.apply(A, b)
if mesh.num_cells() < 16:
print(A.array())
print(b.array())


With access to the elements in A through a numpy array, we can easily perform computations on this matrix, such as computing the eigenvalues (using the eig function in numpy.linalg). We can alternatively dump A.array() and b.array() to file in MATLAB format and invoke MATLAB or Octave to analyze the linear system. Dumping the arrays to MATLAB format is done by

import scipy.io
scipy.io.savemat('Ab.mat', {'A': A.array(), 'b': b.array()})


Writing load Ab.mat in MATLAB or Octave will then make the array variables A and b available for computations.

Matrix processing in Python or MATLAB/Octave is only feasible for small PDE problems since the numpy arrays or matrices in MATLAB file format are dense matrices. FEniCS also has an interface to the eigensolver package SLEPc, which is the preferred tool for computing the eigenvalues of large, sparse matrices of the type encountered in PDE problems (see demo/documented/eigenvalue/python/ in the FEniCS/DOLFIN source code tree for a demo).