$$ \newcommand{\dt}{\Delta t} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\x}{\boldsymbol{x}} \renewcommand{\u}{\boldsymbol{u}} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \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} \newcommand{\renni}[2]{\langle #2, #1 \rangle} $$
Solving PDEs in Minutes -
The FEniCS Tutorial Volume I

 

 

 

Choosing a test problem

The Poisson problem (2.1)--(2.2) has so far featured a general domain \( \Omega \) and general functions \( \ub \) for the boundary conditions and \( f \) for the right-hand side. For our first implementation we will need to make specific choices for \( \Omega \), \( \ub \), and \( f \). It will be wise to construct a problem where we can easily check that the computed solution is correct. Solutions that are lower-order polynomials are primary candidates. Standard finite element function spaces of degree \( r \) will exactly reproduce polynomials of degree \( r \). And piecewise linear elements (\( r=1 \)) are able to exactly reproduce a quadratic polynomial on a uniformly partitioned mesh. This important result can be used to verify our implementation. We just manufacture some quadratic function in 2D as the exact solution, say $$ \begin{equation} \tag{2.12} \uex(x,y) = 1 +x^2 + 2y^2\tp \end{equation} $$ By inserting (2.12) into the Poisson equation (2.1), we find that \( \uex(x,y) \) is a solution if $$ f(x,y) = -6,\quad \ub(x,y)=\uex(x,y)=1 + x^2 + 2y^2,$$ regardless of the shape of the domain as long as \( \uex \) is prescribed along the boundary. We choose here, for simplicity, the domain to be the unit square, $$ \Omega = [0,1]\times [0,1] \tp$$ This simple but very powerful method for constructing test problems is called the method of manufactured solutions: pick a simple expression for the exact solution, plug it into the equation to obtain the right-hand side (source term \( f \)), then solve the equation with this right-hand side and try to reproduce the exact solution.

Tip: Try to verify your code with exact numerical solutions! A common approach to testing the implementation of a numerical method is to compare the numerical solution with an exact analytical solution of the test problem and conclude that the program works if the error is "small enough". Unfortunately, it is impossible to tell if an error of size \( 10^{-5} \) on a \( 20\times 20 \) mesh of linear elements is the expected (in)accuracy of the numerical approximation or if the error also contains the effect of a bug in the code. All we usually know about the numerical error is its asymptotic properties, for instance that it is proportional to \( h^2 \) if \( h \) is the size of a cell in the mesh. Then we can compare the error on meshes with different \( h \)-values to see if the asymptotic behavior is correct. This is a very powerful verification technique and is explained in detail in the section Computing convergence rates. However, if we have a test problem for which we know that there should be no approximation errors, we know that the analytical solution of the PDE problem should be reproduced to machine precision by the program. That is why we emphasize this kind of test problems throughout this tutorial. Typically, elements of degree \( r \) can reproduce polynomials of degree \( r \) exactly, so this is the starting point for constructing a solution without numerical approximation errors.