$$ \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}[2]{\langle #1, #2 \rangle} $$

 

 

 

Extensions: Improving the Poisson solver

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.

Refactoring the Poisson solver

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.

A more general solver function

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.

Writing the solver as a Python module

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.

Verification and unit tests

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

Regarding the last point, 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.

Tip: Print messages in test functions. The 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.

Tip: Debugging with iPython. One can easily enter iPython from a Python script by adding the following line anywhere in the code:

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.

Parameterizing the number of space dimensions

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)