I wanted to verify my FEM-program by applying the method of manufactured solutions, while solving the Poisson equation in two dimensions using the continuous Galerkin method
$$-\nabla^2u=f$$ with $$u=\sin(\pi x)\sin(\pi y)$$ resulting in $$f=2\pi^2\sin(\pi x)\sin(\pi y)$$ My system is defined on a square with $\Omega=[0,1]^2$, and dirichlet boundary conditions: $u_{\partial\Omega}=0$
Then I looked at the convergence rate of the $H^1$-seminorm for elements from second to eighth order at different grid densities, resulting in
+--------------+---------+----------+----------+------------+-------------+--------------+---------------+
| cell numbers | P2 | P3 | P4 | P5 | P6 | P7 | P8 |
+--------------+---------+----------+----------+------------+-------------+--------------+---------------+
| 4 | 11.04 | 339 | 852.46 | 10595.67 | 137926.51 | 6300733.33 | 69993038.13 |
| 16 | 3.97 | 62.5 | 12.8 | 14.29 | 56.35 | 335.63 | 295.60 |
| 64 | 4 | 125.47 | 12.81 | 14.02 | 56.95 | 350.12 | 299.23 |
| 256 | 4 | 251.15 | 12.81 | 13.96 | 57.10 | 303.92 | 12.77 |
| 1024 | 4 | 502.4 | 12.8 | 13.94 | 57.24 | 6.82 | 1.16 |
+--------------+---------+----------+----------+------------+-------------+--------------+---------------+
Some things here I (think I can) explain, others not. On the one hand, the large value in the first line results from the fact that I start from $u_0=0$, while for every other step I reuse the result from the last step, which has been interpolated onto the new grid. Thus I expect the correct convergence rate showing up there, except in the first step.
Furthermore, I assume that the decreasing values for elements of seventh and eighth order show that the interpolation already is quite good for low grid densities. Is that correct?
The things I do not understand here:
Why is the convergence rate for elements of third order so much higher (and increasing for increasing density) than for the elements of fourth and fifth order (partially also for higher orders)? Is that by accident, i.e. third-order elements simply interpolate the test function quite well, or something else?
Why is there a big step in convergence rates from second order to fourth order (skipping third order), but almost no difference between the fourth and fifth order? Again, a step from fifth to sixth order and from sixth to seventh, but almost no difference is visible between the seventh and eighth order.
How can I explain those things? (I hope there is no error in the program itself...)
Edit: The values of the $H^1$-seminorm itself are
+--------------+----------+----------+----------+-----------+-----------+-----------+-----------+
| cell numbers | P2 | P3 | P4 | P5 | P6 | P7 | P8 |
+--------------+----------+----------+----------+-----------+-----------+-----------+-----------+
| 4 | 3.485e-1 | 1.135e-2 | 4.514e-3 | 3.631e-4 | 2.790e-5 | 6.107e-7 | 5.497e-8 |
| 16 | 8.820e-2 | 7.335e-4 | 2.897e-4 | 1.172e-5 | 4.451e-7 | 4.343e-9 | 2.131e-10 |
| 64 | 2.210e-2 | 4.624e-5 | 1.823e-5 | 3.684e-7 | 6.994e-9 | 3.293e-11 | 8.274e-13 |
| 256 | 5.527e-3 | 2.896e-6 | 1.141e-6 | 1.157e-8 | 1.094e-10 | 2.981e-13 | 7.630e-14 |
| 1024 | 1.382e-3 | 1.811e-7 | 7.135e-8 | 3.617e-10 | 1.708e-12 | 1.183e-13 | 1.428e-13 |
+--------------+----------+----------+----------+-----------+-----------+-----------+-----------+
I am using the library deal.II for calculating the values, and use the element FE_Q
as basis element. It switches from lagrange polynomials (up to order 2) with equidistant support points to Gauss-Lobatto-polynomials for higher orders (also refer to https://www.dealii.org/developer/doxygen/deal.II/classFE__Q.html). Could that be the reason for the observed convergence behaviour?