Integrate DG into non-nested MG

[jh66637]

This PR aims to integrate DG into the non-nested MG implementation and is a split off of #15148.

@peterrum @fdrmrc

[jh66637]

@peterrum I think I included all requested changes :)

[peterrum]

First set of comments. Thank you for working on this!

[jh66637]

@fdrmrc This should be ready. Maby you can have a look and test it :)

Regarding an extension to multiple components: Probably only these lines have to be changed
https://github.com/jh66637/dealii/blob/7f5b41616e4a3b21426d05122394bb3e4d638dfb/in[…]ude/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
considering
https://www.dealii.org/developer/doxygen/deal.II/classFiniteElement.html#ae2ea16b60a6fc644a9bc7097703a53e8.

@peterrum Thanks for the help :)

[peterrum]

Nice feature of 9.5 :)

[peterrum]

@jh66637 Thank you for the very nice work! These modifications are also a big step towards supporting multiple components in the non-nested MG implementation needed by @fdrmrc !

Let's wait for the review of another person, since I have been partly involved in the development. Please note there might be still some changes coming depending on the requirements of multi-component systems.

[kronbichler]

This is a complicated algorithm and I would like to ask if I understood the overall idea correctly:

• We are primarily discussing the case of a prolongate operation that goes into a DG space.
• The main problem with a DG space is that several elements might request the evaluation at the same point. Rather than removing duplicates among the points as in the FE_Q case (as implemented via Add two level transfer operator between non-nested levels #15141), we need to take into account all values, as they in fact might be different on different sides for a discontinuous basis. This is done by constructing an auxiliary continuous finite element space, which has the same support points as the underlying DG space, which is underlying the RemotePointEvaluation process.

So in some sense, this could be seen as a performance optimization to reduce the work to the points actually needing the information.

Is this understanding correct? If no, could you please help me understand where in the above chain my misunderstanding is? Use the comments below as guidance where my misunderstanding could come from.

If my understanding is correct, I am confused why we need the complicated code. The point I want to understand is what would happen with FE_DGQArbitraryNodes<dim>(QGauss<1>(degree+1)). As I explain more thoroughly in one of my comments, also the usual FE_DGQ could be interpreted as an arbitrary-nodes interpolation. As the grids will in general not be matching, I see no reason why we would stick to FE_DGQ and not the more "accurate" one with Gauss quadrature points. This question evolved during the ~3 hours when reviewing this PR, spurred by offline discussions with @peterrum and the special cases you need for degree 0: In a sense, the degree 0 case is the prototypical FE_DGQArbitraryNodes(QGauss<1>(degree+1)) element, with unique points on the grid. My question would then be if, for all kinds of FE_DGQ* elements, we should not always interpolate into some arbitrary support points that are rich enough first, and only later switch into the actual basis representation?

And finally, I would be wondering if/how we would want to proceed if we wanted to make some consistent projection from one grid to the other, say, by integration on cut sub-entities of the mesh where each entity has polynomial solutions from both sides as generated by CGAL intersections. Where would one put that in this framework?

I realize my request may sound overly picky and my questions might entirely miss the intent you have, but in fact I highly appreciate the work you've been doing, this is excellent work and helpful in a variety of contexts.

[kronbichler]

Do we need some additional requirements on the two spaces? Like the number of dofs_per_cell on the base element being the same? Just because two FE spaces have the same degree does not mean they are compatible one-to-one. As far as I can see in the code below, you only check the degree of the finite element before you create the dof_handler_sp, without ever looking at the unknowns. What would happen if I used this with FE_DGP, https://dealii.org/developer/doxygen/deal.II/classFE__DGP.html?

[jh66637]

Thank you for pointing out the missing Asserts.
I added Asserts to ensure the function can only be used with FE_DGQ and FE_Q. I also added an Assert to ensure the same n_dofs_per_cell which is indeed needed for now (and has to be deleted later when the function is extended such that it can be used to identify DoF indices for the same support point in the context of multiple components).

7b8ee1b

[kronbichler]

This assumes sp_indices.size() == dof_indices.size() and the same position of support points, right? And that the dof_handler actually have an underlying FiniteElement with support points? I think we need to make additional assertions that the finite element spaces satisfy the underlying assumptions of this algorithm. This is an extension of the comment above, and to test what I mean, I suggest to consider the case of FE_DGQLegendre, https://dealii.org/developer/doxygen/deal.II/classFE__DGQLegendre.html and see if that one also works.

[jh66637]

I will update the asserts tomorrow, I want to take the time to test this as suggested.

[kronbichler]

I do not understand why this is necessary, given the comment of FE_DGQ<dim>(0) you make above. There are also DG elements with unique support points, or in fact they could all be considered as having unique support points. My reasoning goes like this: The DG polynomial space is equivalent to the Lagrange basis in the Gaussian quadrature points, which are all located strictly inside the reference element. The reason why you do the averaging above (which took me a long time to understand, it is complicated) is that you assume the Gauss-Lobatto support points, which is however an arbitrary choice in a DG method. By construction, a DG element does not impose continuitiy in any way, so whether or not two support points coincide can be considered by chance but not systematic. I do not want to dispute that there are cases where achieving continuity in the interpolation between grids makes sense, but I would like to see some theoretic justification. I do really appreciate the effort you put into making the matching of unknowns work, but as I told to @peterrum offline, my question is what you would like to happen if I pass in FE_DGQArbitraryNodes<dim>(QGauss<1>(degree + 1));. I summarize my point in the general introductory comment to this review.

[jh66637]

I completely agree and dependent on which direction we agree to move I will change this 👍

[kronbichler]

As a note, several of the above theoretical questions do not need to be addressed immediately and can be considered part of the more general issue #15148, right now I would be happy with some clarifications to see the overall direction more clearly.

[jh66637]

@kronbichler Thank you for reviewing, I know the PR is quite lengthy and I highly appreciate it. I will try to answer your questions. If something remains unclear, please let me know.

Before answering all the questions you left above, in detail, I would also like to discuss the main picture. This is because I see your point and would like to know in which direction we are targeting.

During the development of the code, I actually only thought about FE_DGQ. Reading through your comments, I think it is valid to see this approach as a performance optimization for FE_DGQ (even though a lot of the code can be reused for multiple components in continuous FEM).

The motivation is as follows: In case of continuous FEM it is possible to remove the duplicate DoF indices and directly access the DoF values in the vector. If we want to directly access the DoF values from a given vector for DG (instead of using additional evaluations) we have to scan for duplicate DoF indices corresponding to the same support point. The latter holds only true if support points are not strictly located inside elements.

Interpolating into arbitrary support points first and switching to the basis representation afterward always requires the evaluation of shape functions instead of directly accessing components of a vector, and using the non-nested version on a nested (globally coarsened grid) would yield different results compared to directly applying the nested version. Apart from this, it has the advantage that it is possible to skip the scanning for indices at support points.

I hope this clarifies the intent. I don't know in which direction we want to move here. One suggestion would be to keep the code as a fast path for the special case with FE_DGQ and take a more general route in other cases, but I am also open to doing it differently.

Regarding where to place consistent projection using CGAL: The code will look quite different, so I don't know where to place it and how interfaces should ideally look like. @peterrum and @fdrmrc: Did you think about the interfaces already?

[peterrum]

@peterrum and @fdrmrc: Did you think about the interfaces already?

I let @fdrmrc comment on this topic, since he is currently working on this. But I think we agreed on an approach where one can configure the two-level transfer operator via an AdditionalData. I believe that although the the point-wise and the l2-approaches differ, they will share significant part of the code: the solution will need to be evaluated at points (support points of the fine mesh for point-wise; quadrature points of the fine mesh; quadrature points on intersections created by CGAL -> #15165).

[kronbichler]

Thank you @jh66637 for the extensive explanation. I am not sure either about the direction to take, I think your point is valid in terms of support points and "optimization" possibilities.

[kronbichler]

Interpolating into arbitrary support points first and switching to the basis representation afterward always requires the evaluation of shape functions instead of directly accessing components of a vector, and using the non-nested version on a nested (globally coarsened grid) would yield different results compared to directly applying the nested version.

I think we would still get the same result, because the prolongation operation that gets implemented is an embedding for nested meshes: No matter what the basis representation on the fine scale is, it will always be able to exactly represent the coarse scale. And in fact, any p+1 points per direction suffice to exactly determine a polynomial of degree p (or two polynomials in the nested case). In the non-nested case, the arbitrary nodes have the advantage of being a bit more accurate in the sense that the represent polynomials of two degrees higher QGauss versus QGaussLobatto, but you pay it with additional evaluations as you observed. And then, the main application of the non-nested transfer is for the case where the solution to be transferred is not smooth, so the higher rate of Gauss versus Gauss-Lobatto is not as important. As you have observed in your applications, it is hard to avoid using consistent integration patches on the intersections to guarantee polynomial representation.

Anyway, I suggest to proceed with the present PR, do the small tidy-up and address possible algorithmic strategies later.

[fdrmrc]

I also think the present version could be used as a fast path for FE_DGQ elements.

I let @fdrmrc comment on this topic, since he is currently working on this.

As written, what we were thinking about was indeed an AdditionalData that allows to configure a transfer operator which can be: interpolation, projection or projection_exact, where the latter computes intersections with CGAL.

I am working on projection (without CGAL) right now and most of the work currently available inside reinit() can be readily used also for classical projection. In essence, the only addition required in this case is the matrix-free setup that can be done at the end of the reinit().

I haven't thought about the projection_exact interface yet, so I can't really comment, but I agree with @jh66637 that the code will look quite different (it will depend on #15165) so it will probably have to be handled separately.