Convergence rate

[Félix Bunel]

Hello everyone.

I'm having some trouble to understand the convergence rate i'm observing in my code.

Here is what i'm solving :

- I'm in 2D on a round mesh.

- I'm solving a simple Poisson equation on this mesh for a variable named Phi the solution is known for this and is 1-x^2-y^2

- With this solution phi I'm then solving a Stokes equation that has special terms that depends on phi.

For the stokes problem i'm using the usual mixed fe element as such :

stokes_fe(FE_Q<2>(2), 2,
          FE_Q<2>(1), 1),

So second order for the speeds and first order for the pressure (just like in the boussinesq problem from the tutorials.

For the poisson/phi problem, i'm using FE_Q<2>  also

initialized as such :

phi_fe (2),

So second order.

For a special case, I have a known solution which is u=v=0 and p of the form 1-x^2-y^2 (just like phi)

And i have solved this on multiple refinement cycle which gives me different number of dofs and cellsize.

The thing is, when I plot the error as a function of the maximum cell diameter, I get a quadratic convergence rate for the speed, and a linear rate for P and phi.

What surprises me is that I don't have a quadratic convergence rate for my poisson/\phi problem even though I used

phi_fe (2),

My norm is the L2 norm which is simply sqrt(sum(error**2))). I'm computing it with python after exporting the solution as a gpl file.

Here are the graphs :

I have also tried

phi_fe (1),

But in this case I don't even have convergence for the speed.

In my integration part, I have always used    

QGauss<2>  quadrature_formula(4);

which is equal to degree+2.

So here are my questions :

1- Is this convergence rate the one i'm expecting for a poisson equation (phi problem) ?

2- If no, any idea why I don't have the correct convergence rate ?

3- In the tutorial step 31, a degree of 2 is used for the temperature (which is very similar to my phi), is there some reason for that choice (just like using degree+1 for the speed and degree for the pressure in a stokes problem...)

Thanks in advance.

[Bruno Blais]

A quick question, since you are working on a sphere, are you specifying a mapping of the same order as your phi?

[Félix Bunel]

Thanks for your answer !

I'm working in 2D on a disk !

This is my mesh on the 4 cycle of refinement (which I never use in practice because it's not refined enough for what I want to do).

I don't think I have a mapping object...

What i do to initalize my system on phi is :

//On génère le dofhandler
    phi_dof_handler.distribute_dofs (phi_fe);

    //On génère les conditions aux bords
    phi_constraints.clear();
    DoFTools::make_hanging_node_constraints(phi_dof_handler, phi_constraints);
    VectorTools::interpolate_boundary_values(phi_dof_handler,
                                             0,
                                             ZeroFunction<2>(1),
                                             phi_constraints
                                            );

    //On créer le sparsity pattern
    DynamicSparsityPattern dsp(phi_dof_handler.n_dofs());
    DoFTools::make_sparsity_pattern (phi_dof_handler, dsp);
    phi_sparsity_pattern.copy_from(dsp);

    //On resize les objets à la bonne taille
    phi_system_matrix.reinit (phi_sparsity_pattern);
    phi_solution.reinit (phi_dof_handler.n_dofs());
    phi_old_solution.reinit (phi_dof_handler.n_dofs());
    phi_system_rhs.reinit (phi_dof_handler.n_dofs());

    //On affiche les infos intéressantes
    std::cout << YELLOW << "     Nombre de degrés de liberté pour les angles : "
              << RESET  << WHITE  << phi_dof_handler.n_dofs() << RESET
              << std::endl;

[luca.heltai]

Dear Felix,

by any chance, did you take a look at step 10?

https://www.dealii.org/current/doxygen/deal.II/step_10.html

This step explains a little bit what to do when you want to solve on curved domains with high order finite elements. In particular, you need to ensure that the mapping you are using (from the reference element to the current element) is at least of the same order of your finite element method, if you want to achieve the correct convergence rates. In your case, you are using a MappingQ1 (the default mapping, if you don’t specify anything). With that, you won’t go beyond 2nd order in L2. If you use MappingQ2, you should see again the optimal convergence rates.

Best,
Luca.

> On 10 Oct 2019, at 14:11, Félix Bunel <bunel...@gmail.com> wrote:
>
> Thanks for your answer !
>
> I'm working in 2D on a disk !
>
> This is my mesh on the 4 cycle of refinement (which I never use in practice because it's not refined enough for what I want to do).
>

> <Auto Generated Inline Image 1.png>

> The thing is, when I plot the error as a function of the maximum cell diameter, I get a quadratic convergence rate for the speed, and a linear rate for P and phi.
> What surprises me is that I don't have a quadratic convergence rate for my poisson/\phi problem even though I used
> phi_fe (2),
>
> My norm is the L2 norm which is simply sqrt(sum(error**2))). I'm computing it with python after exporting the solution as a gpl file.
> Here are the graphs :
>
>
>

> I have also tried
> phi_fe (1),
>
> But in this case I don't even have convergence for the speed.
>
>
>

> In my integration part, I have always used
> QGauss<2> quadrature_formula(4);
> which is equal to degree+2.
>
> So here are my questions :
>
> 1- Is this convergence rate the one i'm expecting for a poisson equation (phi problem) ?
>
> 2- If no, any idea why I don't have the correct convergence rate ?
>
> 3- In the tutorial step 31, a degree of 2 is used for the temperature (which is very similar to my phi), is there some reason for that choice (just like using degree+1 for the speed and degree for the pressure in a stokes problem...)
>
> Thanks in advance.
>

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/966a3637-8ca8-462a-bf7a-621746c9c6ad%40googlegroups.com.
> <Auto Generated Inline Image 1.png>

[Félix Bunel]

Hi thank you very much !

That was the problem I believe ! THANKS A LOT ! Now the error go down quadratically at the beginning.

Then it goes up again linearly which i find very strange... I believe it's because I haven't set enough precision in the solving process or maybe i'm reaching machine precision...

> To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.

[Félix Bunel]

I found why it's not going further down ! It was just that the GPL file saved is with a precision of 6 digits which corresponds with the error per dofs when I norm my error plot by the solution.

Thanks for the help everyone !