3
$\begingroup$

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?

Improve this question
$\endgroup$
1
  • $\begingroup$ Could you be more specific about what the values are in both of these tables? It's really unclear to me. Is the first table a ratio between consecutive error norms, possibly?? Is the second table the seminorm of the solution itself, or of the error? Formulas for these would be very helpful. $\endgroup$
    – LedHead
    Jan 22, 2019 at 0:10

3 Answers 3

Reset to default
4
$\begingroup$

There is a bug in your code :-) First, going down in each column, you have cases where the error becomes larger again, and this should clearly not happen because the finite element spaces are nested going from top to bottom. Secondly, moving left to right the error should also decrease because the finite element spaces are nested, but again this is not happening. Finally, in most columns the error does not appear to converge to zero but to some finite value -- this should also clearly not happen.

So if the table shows errors, then there is a bug somewhere. If, in contrast, the table shows the $H^1$ seminorm of the solution rather than the error, then there is still a bug somewhere because then going top to bottom and left to right should still result in convergence towards the exact value, but your values are all over the place.

Improve this answer
$\endgroup$
4
  • $\begingroup$ The table shows the convergence rate of the $H^1$-seminorm, i.e. the convergence rate from 4 cells to 16 cells in the case of P2-elements is ~4, while for P8-elements it is 295. If required, I also can add the value for the $H^1$-seminorm itself. $\endgroup$
    – arc_lupus
    Jan 21, 2019 at 18:14
  • $\begingroup$ I added the values from the $H^1$-seminorm itself to the question $\endgroup$
    – arc_lupus
    Jan 21, 2019 at 19:20
  • $\begingroup$ I am using FE_Q-elements from deal.II for the calculation. Could there be a relation between the behaviour and the elements? $\endgroup$
    – arc_lupus
    Jan 22, 2019 at 10:20
  • $\begingroup$ I think I'm not quite clear about what you are showing. You say it's the H1 seminorm, but the seminorm of what? The error? The solution? $\endgroup$
    – Wolfgang Bangerth
    Jan 22, 2019 at 20:43
1
$\begingroup$

How are you actually computing error norm and convergence rate ?

If I use your error table, I get the expected result. E.g.

P2, 256 and 1024 cells

In [1]: log(5.527e-3/1.382e-3)/log(2)                                           
Out[1]: 1.9997389968259023

P4, 256 and 1024 cells

In [4]: log(1.141e-6/7.135e-8)/log(2)                                           
Out[4]: 3.9992415517002806

There is something weird going on with your P3 results. You must double check your code.

Update: I have a deal.II code that solves almost the problem here

https://github.com/cpraveen/fembook/tree/master/deal.II/ex04

Modifying this for your particular problem with degree = 3, I get

cells dofs       L2             H1       
    4   49 1.359e-03    - 2.668e-02    - 
   16  169 8.812e-05 3.95 3.376e-03 2.98 
   64  625 5.564e-06 3.99 4.233e-04 3.00 
  256 2401 3.486e-07 4.00 5.295e-05 3.00 
 1024 9409 2.180e-08 4.00 6.620e-06 3.00
Improve this answer
$\endgroup$
3
  • $\begingroup$ Do you have any suggestions where to start? $\endgroup$
    – arc_lupus
    Jan 22, 2019 at 9:37
  • 1
    $\begingroup$ @arc_lupus Difficult to say, I have trouble finding my own bugs :-( I added link to a deal.II code that shows no problem for P3. Maybe it helps you to debug your code. $\endgroup$
    – cfdlab
    Jan 22, 2019 at 13:57
  • $\begingroup$ I found the problem, thanks to your code... Thanks! $\endgroup$
    – arc_lupus
    Jan 22, 2019 at 15:08
0
$\begingroup$

I found the explanation for that behaviour: For the calculation of the integral over the differences I used an element of third order, which was fixed. This explains the high convergence rate for P3, and the lower rate for all other elements. When replacing the fixed order with a more suitable order (such as FE_degree * 2 + 1), my convergence rate becomes

P2    4
P3    8
P4    16
P5    32
P6    64
P7    128
P8    256

as expected.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Ah, that's better :-) You could also have used the VectorTools::integrate_difference() function if you happen to know the exact solution. $\endgroup$
    – Wolfgang Bangerth
    Jan 22, 2019 at 20:44
  • $\begingroup$ According to https://www.dealii.org/developer/doxygen/deal.II/step_7.html I still have to use a quadrature norm when calling VectorTools::integrate_difference() with a certain degree. In the example a fixed degree was used, while cpraveen used 2*fe.degree + 1. I strongly assume that that was the problem. $\endgroup$
    – arc_lupus
    Jan 23, 2019 at 10:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.