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

- has a name starting with
`test_`

, - has no arguments, and
- implements a test expressed as
`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)
```