Error in the second derivative of solution...?

[hank...@gmail.com]

Dear All,

I'm trying to calculate some problem using deal.II.

But I have some problem with it....

First, what I am calculating is q value as you can see above(that includes the second derivatives of solution that is psi)...

psi is the solution of Laplace equation with coefficient.

In any case, I calculated q with FE_Q<dim> and attached the several results with various order of FE_Q.

However, as you can see in the picture I attach, even though the shape of profile is almost same, the results(q values itself) are very different depending on the order of FE_Q.

For instance When I use 15314 experimental data for the input of my deal.II code, q has minus sign with Q3 or Q4 finite element and plus sign with Q5 or Q6 element.

However, according to the theory of plasma physics, q can't be minus value unlike q values of 15314_Q3, Q4 and can't be so high unlike q values of 15314_Q6 

So, I thought Q5 was best choice for the calculation of q.

However, When I input 13256 experimental data for the calculation, it seems that Q6 element is best.

And then I tried to input many other experimental data, and I can't determine yet which order to use....

I guess this problem comes from the second derivatives of solution.

So my question is....

Is there any Finite element that is most appropriate for the second derivative?

Thank you 

Kyusik.

[Wolfgang Bangerth]

On 05/14/2017 01:57 AM, hank...@gmail.com wrote:
>
> I guess this problem comes from the second derivatives of solution.

I don't know whether the problem really comes from that, but you can't take
the second derivative of a finite element solution without problem. The issue
is that it exists, of course, in the interior of elements. But it is a
delta-function like thing at the interfaces between cells because the finite
element solution is continuous but not continuously differentiable across the
interface (it has a "kink"), and consequently the second derivative consists
of delta functions. I suppose you are not taking them into account in your
formulation.

If you do need second derivatives, you need to either use a C^1 element (quite
difficult to do), or use a mixed formulation in which you rewrite the equation
-Delta q = f into
u + nabla q = 0
-div u = -f
and then use a continuous element for the variable u. This way u=-grad q,
and Delta q = -div u, which you can compute because u is continuous.

Best
W.

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

[hank...@gmail.com]

Dr. Wolfgang Bangerth,

Thank you very much for your kind reply.

But, could I have one more question....?

I don't know whether the problem really comes from that, but you can't take 
the second derivative of a finite element solution without problem. The issue 
is that it exists, of course, in the interior of elements.

I just had a good idea for the calculation q without second derivative yesterday.

But I need the surface integral of divergence of f(R,solution)*grad_solution

where f is function of position and solution.

So my new question is , in deal.II, could I get local divergence of certain function to get surface integral?

Thank you.

Kyusik.

[Jean-Paul Pelteret]

Dear Kyusik,

So my new question is , in deal.II, could I get local divergence of certain function to get surface integral?

Yes it is. You can extract the divergence of a function using, for example, the FEValuesView function get_function_divergences. 

However, you'll note that this function only makes sense (and its therefore only defined) for a vector field. Your solution field is scalar though, which presents a bit of an issue. On the face of it you'd have to first project the gradient of the solution (scaled by the function "f") onto an axillary vector-valued FE space and then compute the divergence of that projected field. 

I think that what you could alternatively do is to expand your integrand into a gradient term and a Laplacian term. You can quite easily point-wise compute these individual quantities using the scalar FEValueViews functions get_function_gradients and get_function_laplacians, or those equivalently defined for FEValuesBase.

Does this help?

Regards,

Jean-Paul

[hank...@gmail.com]

Dear Jean-Paul Pelteret,

Thank you very much for your help. It is very helpful.

However, you'll note that this function only makes sense (and its therefore only defined) for a vector field. Your solution field is scalar though, which presents a bit of an issue. On the face of it you'd have to first project the gradient of the solution (scaled by the function "f") onto an axillary vector-valued FE space and then compute the divergence of that projected field. 

I'm sorry to tell you that actually I don't know how to project grad_solution into the vector-valued FE space....

Do you happen to know some example tutorial code that includes it....?

Thank you.

Kyusik.   

[Jean-Paul Pelteret]

Dear Kyusik,

 

I'm sorry to tell you that actually I don't know how to project grad_solution into the vector-valued FE space....

Do you happen to know some example tutorial code that includes it....?

Off the top of my head, no I don't know of an example tutorial that demonstrates it. But you're effectively solving a least-squares problem for the solution gradient so its relatively straight-forward to do:

1. Create a vector-valued finite element system (maybe an FESystem<dim>(FE_Q<dim>(1),dim) ) and use it to distribute DoFs onto a new DoFHandler.

2. Create a mass matrix "M" with the new DoFHandler.

3. Create a RHS vector "f" of the form f_i = \int_\Omega f(x) \phi_i(x) dx . Here f(x) = \grad \psi (x) is the gradient of the solution evaluated at point x and \phi_i is a vector-valued shape function. So you need to obviously compute the gradient of the solution at all quadrature points when doing integration - you would do this with the function I mentioned previously.

4. Solve the linear system M.d = f, where "d" represents the (continuous) gradient of the solution to the scalar problem.

There's some useful information in the introduction to the VectorTools namespace (namely on projection and the form of the RHS vector), so you can look there as well for more information. 

As I said before, you could avoid these complexities by expanding your original equation into the two separate terms. However, the benefit to the projection approach is that you can evaluate this term for one order lower polynomial since you don't need to directly compute second derivatives. With the projection one might incur some loss of accuracy of the gradient field, but I'm not quite sure to what extent this would affect your calculation. 

Regards,

Jean-Paul

[hank...@gmail.com]

Dear Jean-Paul,

Thank you again for a detailed explanation.

However, please excuse my ignorance... 

1. Create a vector-valued finite element system (maybe an FESystem<dim>(FE_Q<dim>(1),dim) ) and use it to distribute DoFs onto a new DoFHandler.

2. Create a mass matrix "M" with the new DoFHandler.

3. Create a RHS vector "f" of the form f_i = \int_\Omega f(x) \phi_i(x) dx . Here f(x) = \grad \psi (x) is the gradient of the solution evaluated at point x and \phi_i is a vector-valued shape function. So you need to obviously compute the gradient of the solution at all quadrature points when doing integration - you would do this with the function I mentioned previously.

What makes me confused here is that, according to your explanation, two different DoFHandler could be used simultaneously.

I mean let's say dof_handler1 is for solution and dof_handler2 is for M

To get \grad\psi, I have to use dof_handler1, however, to assemble matrix M and vector f I have to use dof_handler2...

So, in the assemble_system() I don't know whether I have to use 

cell= dof_handler1.begin_active(), endc = dof_handler1.end();

or 

cell= dof_handler2.begin_active(), endc = dof_handler2.end()

Also, It seems that I have to use two different FEValues<dim> one for solution the other for the vector...

In the assemble_system() for M and f, can I use these two things simultaneously...?

Thank you

Kyusik.

 

[Jean-Paul Pelteret]

Hi Kyusik,

What makes me confused here is that, according to your explanation, two different DoFHandler could be used simultaneously.

Yes, that it correct. You can convert between cell iterators. So I do the following:

          typename Triangulation<dim>::active_cell_iterator

          cell = triangulation.begin_active(),

          endc = triangulation.end();

          for (; cell!=endc; ++cell)

          {

            const typename DoFHandler<dim>::active_cell_iterator cell_1 (&triangulation, cell->level(), cell->index(), &dof_handler_1);

            typename DoFHandler<dim>::active_cell_iterator cell_2 (&triangulation, cell->level(), cell->index(), &dof_handler_2);
          ...}

and then follow the usual process of FEValues reinitialisation

 

Also, It seems that I have to use two different FEValues<dim> one for solution the other for the vector...

Yes you would need to. The one relates to the scalar FE system and the other to the new vector FE system.

 

In the assemble_system() for M and f, can I use these two things simultaneously...?

Sure. You use one for assembly and one for the extraction of gradient data from the scalar solution field.

Regards,

Jean-Paul

 

Thank you

Kyusik.

 

[hank...@gmail.com]

 Thank you very much!!

It seems that I can solve my problem.

Thank you!

Kyusik. 

[Jean-Paul Pelteret]

Thats great! I'm glad that you've managed to work it out :-)

Best,

Jean-Paul