hp convergence for non-linear solver.

[Jaekwang Kim]

Hi, all. 

I have asked whether numerical errors also behave ~O(h^p) , even though governing equation is  non-linear and thus, a solver includes iterative method. 

I recall that Dr. Bangerth replied me it is "yes" if we imply small enough tolerance on iteration tolerance. 

I am trying to confirm this testing the code with Picard Iteration, ( Dr. Bangerth's video lecture from 31.5~31.7)

To bring the governing equation again 

Since I want to test the code really has correct h-p convergence, I used the method of manufactured solution. 

To be specific, I assumed solution as 

" u(x,y)="

and this gives me a right hand side function, f ,

f=where z=x^2+y^2

After that I imposed boundary condition as exact with u(x,y) = =g 

Lastly, I compared my numerical solution with manufactured exact solution using Vectortools::integrate_difference. l2norm

template <int dim>

void msurface<dim>::evaluate_error()

{

    Vector<float> difference_per_cell (triangulation.n_active_cells());

    

    VectorTools::integrate_difference (dof_handler,solution,Solution<dim>(),

                                       difference_per_cell,QGauss<dim>(degree+2),VectorTools::L2_norm);

    

    const double L2_error = difference_per_cell.l2_norm();

    std::cout << "   L2_error : " << L2_error << std::endl;

   

    error=L2_error;

}

However, I the error graph I get was as follow...

I think that I didn't make any mistake on my code lines because at least my numerical solution is converging for the case p=1. 

but I couldn't get such convergence for the case p=2, and 3. I have tried to change tolerance number in iterative scheme, but it didn't work. 

Are there other things that I failed to consider? 

Or did I make wrong way to test the code? 

Always thank you for all smart guys here!! 

Jaekwang Kim  

[Bruno Turcksin]

Jaekwang,

On Tuesday, November 1, 2016 at 5:43:36 PM UTC-4, Jaekwang Kim wrote:

I think that I didn't make any mistake on my code lines because at least my numerical solution is converging for the case p=1. 

but I couldn't get such convergence for the case p=2, and 3. I have tried to change tolerance number in iterative scheme, but it didn't work. 

Are there other things that I failed to consider? 

Or did I make wrong way to test the code?

What you did looks good. You should check that you use the right quadrature everywhere. If somewhere in your code you use a low order quadrature, it will work for p=1 but you will have a larger error when you increase the quadrature. You should also look at the solution for p=2 and p=3. This might help you understand what is going on.

Best,

Bruno

[Jaekwang Kim]

Thanks for replying! 

Here I attached solution for p=2, and it seems that it has right share. 

When you said 'right quadrature' dose that mean....quadrature for integration when we assemble system? 

To attach some of my code. 

template <int dim>

void msurface<dim>::assemble_system ()

{

  const MappingQ<dim> mapping (degree);

  const QGauss<dim>  quadrature_formula(degree+2);

  FEValues<dim> fe_values (mapping, fe, quadrature_formula,

                           update_values    |  update_gradients |

                           update_quadrature_points  |  update_JxW_values);

....

I am pretty sure that I am using right quadrature here.... 

but I am not sure whether I am using right quadrature when I compare my numerical solution to exact solution.. 

    VectorTools::integrate_difference (dof_handler,solution,Solution<dim>(),

                                       difference_per_cell,QGauss<dim>(degree+2),VectorTools::L2_norm);

Thanks again for fast reply!! 

2016년 11월 1일 화요일 오후 5시 4분 40초 UTC-5, Bruno Turcksin 님의 말:

[Daniel Arndt]

Jaekwang, 

[...]

    VectorTools::integrate_difference (dof_handler,solution,Solution<dim>(),

                                       difference_per_cell,QGauss<dim>(degree+2),VectorTools::L2_norm);

somegthing that immediately comes to mind is that you don't use a non-default Mapping in here although your boundary is curved.

Make sure that you use an appropriate Mapping everywhere!

Best,

Daniel

[Jaekwang Kim]

Thank for the reply! 

I just fixed my code according to your suggestion. 

Now, I am evaluating my error with following code lines. 

template <int dim>

void msurface<dim>::evaluate_error()

{

    Vector<float> difference_per_cell (triangulation.n_active_cells());

    

    const MappingQ<dim> mapping (degree);

    

    VectorTools::integrate_difference (mapping, dof_handler,solution,Solution<dim>(),

                                       difference_per_cell,QGauss<dim>(degree+2),VectorTools::L2_norm);

    

    const double L2_error = difference_per_cell.l2_norm();

    std::cout << "   L2_error : " << L2_error << std::endl;

   

    error=L2_error;

}

But still. I have following results......

do you have any other idea? 

Thanks again for your time

Regards, 

Jaekwang Kim 

2016년 11월 1일 화요일 오후 5시 19분 35초 UTC-5, Daniel Arndt 님의 말:

[Wolfgang Bangerth]

On 11/01/2016 11:06 PM, Jaekwang Kim wrote:
>
> do you have any other idea?
>

Output the errors per cell and visualize. Is there a pattern? Is the error
localized on the boundary, for example?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Jaekwang Kim]

Thank you for replying me! 

I'd like check the visualized error per cell.

Is there any module in deal.ii that enables this...? 

2016년 11월 2일 수요일 오전 7시 27분 56초 UTC-5, Wolfgang Bangerth 님의 말:

[Wolfgang Bangerth]

On 11/02/2016 07:17 PM, Jaekwang Kim wrote:
>
> I'd like check the visualized error per cell.
> Is there any module in deal.ii that enables this...?
>

Yes, you just need to create a vector with as many entries as there are cell
(e.g., the output of VectorTools::integrate_difference) and attach that to the
DataOut object you have. step-27 shows how to do that, but it in essence works
just like for vectors that have as many entries as there are DoFs.

[Jaekwang Kim]

Thanks for replying me ! 

I moved to squared domain to avoid other possible error when I imply boundary condition, but I still have same manufactured solution e^(-x^2-y^2)

I have to say that I have some changed error curve after moving to new mesh domain. 

Which is slightly better, but still unsatisfying results. 

To see error per cell, I referred step 27, and make estimated error data out. 

I used following code lines....

 Vector<float> estimated_error_per_cell (triangulation.n_active_cells());

    

    KellyErrorEstimator<dim>::estimate (dof_handler,

                                        QGauss<dim-1>(3),

                                        typename FunctionMap<dim>::type(),

                                        solution,

                                        estimated_error_per_cell);

In fact, I don't know how kelly error estimator works...

but just checked that it usually works well with Poisson Equations.

Seeing pictures above,,, when p=2 solution looks seemingly good, but estimated error_per_cell should be wrong. 

Since, exact solution the closer points to x=(0,0) (or left down corner in pic) has  sharp gradient, thus the larger error is expected. 

Think estimated error case for p=1 is correct. 

Now, I confirmed that the code is not working at higher degree. (p>=2)...

but it is not because of boundary condition...since I have four boundary (top,bottom,right,left...) but the points which show highest error is some what irrealted to boundary...

Do you have any idea about this? 

Always thank you for your time and advice!

Best regards,

Jaekwang Kim 

2016년 11월 2일 수요일 오후 8시 30분 19초 UTC-5, Wolfgang Bangerth 님의 말: