The FEniCS programs we have written so far have been designed as flat Python scripts. This works well for solving simple demo problems. However, when you build a solver for an advanced application, you will quickly find the need for more structured programming. In particular, you may want to reuse your solver to solve a large number of problems where you vary the boundary conditions, the domain, and coefficients such as material parameters. In this chapter, we will see how to write general solver functions to improve the usability of FEniCS programs. We will also discuss how to utilize iterative solvers with preconditioners for solving linear systems, how to compute derived quantities, such as, e.g., the flux on a part of the boundary, and how to compute errors and convergence rates.
Most programs discussed in this book are "flat"; that is, they are not organized into logical, reusable units in terms of Python functions. Such flat programs are useful for quickly testing ideas and sketching solution algorithms, but are not well suited for serious problem solving. We shall therefore look at how to refactor the Poisson solver from the chapter Fundamentals: Solving the Poisson equation. For a start, this means splitting the code into functions. But refactoring is not just a reordering of existing statements. During refactoring, we also try to make the functions we create as reusable as possible in other contexts. We will also encapsulate statements specific to a certain problem into (non-reusable) functions. Being able to distinguish reusable code from specialized code is a key issue when refactoring code, and this ability depends on a good mathematical understanding of the problem at hand (what is general, what is special?). In a flat program, general and specialized code (and mathematics) are often mixed together, which tends to give a blurred understanding of the problem at hand.
We consider the flat program
ft01_poisson.py
for solving the Poisson problem developed
in the chapter Fundamentals: Solving the Poisson equation.
Some of the code in this program
is needed to solve any Poisson problem \( -\nabla^2 u=f \) on \( [0,1]\times
[0,1] \) with \( u=\ub \) on the boundary, while other statements arise from
our simple test problem. Let us collect the general, reusable code in
a function called solver. Our special test problem will then just be
an application of our solver with some additional statements. We limit
the solver function to just compute the numerical
solution. Plotting and comparing the solution with the exact solution
are considered to be problem-specific activities to be performed
elsewhere.
We parameterize solver by \( f \), \( \ub \), and the resolution of the
mesh. Since it is so trivial to use higher-order finite element
functions by changing the third argument to FunctionSpace, we
also add the polynomial degree of the finite element function space
as an argument to solver.
from fenics import *
import numpy as np
def solver(f, u_D, Nx, Ny, degree=1):
"""
Solve -Laplace(u) = f on [0,1] x [0,1] with 2*Nx*Ny Lagrange
elements of specified degree and u = u_D (Expresssion) on
the boundary.
"""
# Create mesh and define function space
mesh = UnitSquareMesh(Nx, Ny)
V = FunctionSpace(mesh, 'P', degree)
# Define boundary condition
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bc)
return u
The remaining tasks of our initial program, such as calling the solver
function with problem-specific parameters and plotting,
can be placed in a separate function. Here we choose to put this code
in a function named run_solver:
def run_solver():
"Run solver to compute and post-process solution"
# Set up problem parameters and call solver
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
f = Constant(-6.0)
u = solver(f, u_D, 8, 8, 1)
# Plot solution and mesh
plot(u)
plot(u.function_space().mesh())
# Save solution to file in VTK format
vtkfile = File('poisson_solver/solution.pvd')
vtkfile << u
The solution can now be computed, plotted, and saved to file by
simply calling the run_solver function.
The refactored code is placed in a file
ft12_poisson_solver.py.
We should make sure that such a file can be imported (and hence
reused) in other programs. This means that all statements in the main
program that are not inside functions should appear within a test
if __name__ == '__main__':. This test is true if the file is executed as
a program, but false if the file is imported. If we want to run this
file in the same way as we can run ft01_poisson.py, the
main program is simply a call to run_solver followed by a call to
interactive to hold the plot:
if __name__ == '__main__':
run_solver()
interactive()
This complete program can be found in the file ft12_poisson_solver.py.
The remaining part of our first program is to compare the numerical
and the exact solutions. Every time we edit the code we must rerun the
test and examine that error_max is sufficiently small so we know
that the code still works. To this end, we shall adopt unit testing,
meaning that we create a mathematical test and corresponding software
that can run all our tests automatically and check that all tests
pass. Python has several tools for unit testing. Two very popular
ones are pytest and nose. These are almost identical and very easy
to use. More classical unit testing with test classes is offered by
the built-in module unittest, but here we are going to use pytest
(or nose) since that will result in shorter and clearer code.
Mathematically, our unit test is that the finite element solution of
our problem when \( f=-6 \) equals the exact solution \( u=\ub=1+x^2+2y^2 \)
at the vertices of the mesh.
We have already created a code that finds the error at the vertices for
our numerical solution. Because of rounding errors, we cannot demand this
error to be zero, but we have to use a tolerance, which
depends on the number of elements and the degrees of the polynomials
in the finite element basis. If we want to test that the
solver function works for meshes up to \( 2\times(20\times 20) \)
elements and cubic Lagrange elements, \( 10^{-10} \) is an appropriate
tolerance for testing that the maximum error vanishes.
To make our test case work together with pytest and nose, we have to make a couple of small adjustments to our program. The simple rule is that each test must be placed in a function that
test_,assert success, msg.success is a boolean expression that is
False if the test fails, and in that case the string msg is
written to the screen. When the test fails, assert raises an
AssertionError exception in Python, and otherwise runs
silently. The msg string is optional, so assert success is the
minimal test. In our case, we will write assert error_max < tol,
where tol is the tolerance mentioned above.
A proper test function for implementing this unit test in the pytest or nose testing frameworks has the following form. Note that we perform the test for different mesh resolutions and degrees of finite elements.
def test_solver():
"Test solver by reproducing u = 1 + x^2 + 2y^2"
# Set up parameters for testing
tol = 1E-10
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
f = Constant(-6.0)
# Iterate over mesh sizes and degrees
for Nx, Ny in [(3, 3), (3, 5), (5, 3), (20, 20)]:
for degree in 1, 2, 3:
print('Solving on a 2 x (%d x %d) mesh with P%d elements.'
% (Nx, Ny, degree))
# Compute solution
u = solver(f, u_D, Nx, Ny, degree)
# Extract the mesh
mesh = u.function_space().mesh()
# Compute maximum error at vertices
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
error_max = np.max(np.abs(vertex_values_u_D - \
vertex_values_u))
# Check maximum error
msg = 'error_max = %g' % error_max
assert error_max < tol, msg
To run the test, we type the following command:
Terminal> py.test ft12_poisson_solver.py
This will run all functions named test_* (currently only the
test_solver function) found in the file and report the results.
For more verbose output, add the flags -s -v.
We shall make it a habit to encapsulate numerical test problems in unit tests as above, and we strongly encourage the reader to create similar unit tests whenever a FEniCS solver is implemented.
assert statement runs silently when the test passes so users may
become uncertain if all the statements in a test function are really
executed. A psychological help is to print out something before assert
(as we do in the example above) such that it is clear that the
test really takes place.
Note that py.test needs the -s option to show printout
from the test functions.
from IPython import embed; embed()
This line starts an interactive Python session which lets you print and plot variables, which can be very helpful for debugging.
FEniCS makes it is easy to write a unified simulation code that can
operate in 1D, 2D, and 3D. As an appetizer, go back to the
previous programs
ft01_poisson.py
or
ft12_poisson_solver.py
and change the mesh construction from UnitSquareMesh(8, 8) to
UnitCubeMesh(8, 8, 8). Now the domain is the unit cube
partitioned into \( 8\times 8\times 8 \) boxes, and each box
is divided into six tetrahedron-shaped finite elements for
computations. Run the program and observe that we can solve a 3D
problem without any other modifications! (In 1D, expressions must be
modified to not depend on x[1].) The visualization allows you to
rotate the cube and observe the function values as colors on the
boundary.
If we want to parameterize the creation of unit interval, unit square, or unit cube over dimension, we can do so by encapsulating this part of the code in a function. Given a list or tuple specifying the division into cells in the spatial coordinates, the following function returns the mesh for a \( d \)-dimensional cube:
def UnitHyperCube(divisions):
mesh_classes = [UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh]
d = len(divisions)
mesh = mesh_classes[d - 1](*divisions)
return mesh
The construction mesh_class[d - 1] will pick the right name of the
object used to define the domain and generate the mesh. Moreover, the
argument *divisions sends all the components of the list divisions
as separate arguments to the constructor of the mesh construction
class picked out by mesh_class[d - 1]. For example, in a 2D problem
where divisions has two elements, the statement
mesh = mesh_classes[d - 1](*divisions)
is equivalent to
mesh = UnitSquareMesh(divisions[0], divisions[1])
The solver function from
ft12_poisson_solver.py
may be modified to solve \( d \)-dimensional problems by replacing the
Nx and Ny parameters by divisions, and calling the function
UnitHyperCube to create the mesh. Note that UnitHyperCube is a
function and not a class, but we have named it using so-called
CamelCase notation to make it look like a class:
mesh = UnitHyperCube(divisions)