order of cells in two MatrixFree objects on the same Triangulation and Quadrature but different FE spaces

[Denis Davydov]

Dear all,

I am playing around with something in the following setup:

1. I have a single p::d::Tria

2. two different DoFHandler's using scalar FE_Q<dim> of different degrees

3. The same quadrature rule (QGauss) used for both cases

4. Two different MatrixFree objects

I am not 100% sure, but it looks like the iteration over SIMD chunks of cells via two different matrix-free

objects actually leads to the different sequence of cells. I see that when I try to calculate some integral quantities

of my quadrature point data. Namely, i see the difference in the integral

my_integral += quadrature_data(cell,q) * fe_eval.JxW(q);

when I use two different MatrixFree objects to initialize FEEvaluation (obviously with the same n_q_points_1d but different fe_degree).

The order of difference is far beyond any numerical roundoff error.

Is this something expected but not described, AFAICT, in documentation?

Regards,

Denis.

[Denis Davydov]

p.s. the mapping is the same for both cases.

[Martin Kronbichler]

Dear Denis,

p.s. the mapping is the same for both cases.

On Monday, January 29, 2018 at 4:16:48 PM UTC+1, Denis Davydov wrote:

Dear all,

I am playing around with something in the following setup:

1. I have a single p::d::Tria

2. two different DoFHandler's using scalar FE_Q<dim> of different degrees

3. The same quadrature rule (QGauss) used for both cases

4. Two different MatrixFree objects

I am not 100% sure, but it looks like the iteration over SIMD chunks of cells via two different matrix-free

objects actually leads to the different sequence of cells. I see that when I try to calculate some integral quantities

of my quadrature point data.

I guess that's possible. Once there are different DoFHandler objects, the algorithm that collects cells and puts them into batches might change. The documentation does state the following (in FEEvaluation):
https://www.dealii.org/developer/doxygen/deal.II/classFEEvaluation.html
section on "Handling several integration tasks and data storage in quadrature points" towards the end of the general class documentation. It says "For the combination of different fields, including different solution vectors that come from different time steps, it is mandatory that all FEEvaluation objects share the same MatrixFree object."

Did you expect anything in addition?

Best,
Martin

[Denis Davydov]

Hi Martin,

Thanks for confirming that it's indeed possible.

I guess i just missed this part. I always forget that FEEvaluation has much more info than MatrixFree class.

I suppose the right approach here is to use `reinit()` which takes a vector of constraints and DoFHandlers and use it 

in combination with different operators. I guess this might be the only thing to mention in that paragraph.

 

Regards,

Denis

[Denis Davydov]

p.s. I can confirm that using reinit() which takes vectors of DoFHandlers and Constraints solved the problem.

@Martin: nevermind about documentation, there is also a paragraph above which discuss similar things. 

Maybe eventually we need to have another matrix-free tutorial to demonstrate all those options.

Regards,

Denis.