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