Rethink design of OperatorBase regarding different FE types (CG, DG, Hdiv)

[nfehn]

Robin boundary conditions require to evaluate the homogeneous action of a PDE operator on boundary faces also for continuous Galerkin (CG) discretizations, which is new functionality compared to the current code version (where we currently only evaluate inhomogeneous Neumann BCs for CG discretizations).

@richardschu We discussed this topic briefly in person. I create this issue so that we can exchange ideas here (or that others can share their opinion as well).

The question is whether we should rethink the design of OperatorBase regarding the different element types and loop types somewhat basically:

• Which boolean variables do we need / are best to describe the different possible configurations? we currently have evaluate_face_integrals(), is_dg (note that this also includes the Hdiv-case!), (and use_cell_based_face_loops)
• Where should these if-statements be written (on the level of matrix_free->loop/cell_loop() or inside the cell/face/boundary-face integrals)? @kronbichler is it a measureable performance penalty when using matrix_free->loop() but leaving e.g. the function with the face integrals empty?

[richardschu]

From my perspective evaluate_face_integrals() is used to skip the face integration or evaluation stuff completely (probably for performance reasons? @kronbichler ), since it just checks if any evaluation_flags or integration_flags are set.
So I would keep evaluate_face_integrals() as is to keep skipping an otherwise needed empty function for face loop integration.

I think the member variable is_dg in OperatorBase might just be poorly named, could that be? That is_dg might just be renamed to has_interior_face_integrals to also cover the Raviart-Thomas case (or any other non-conforming space that would then result in some interior face integrals being added in some weak form). I think then it would be more intuitive looking at possible extensions in the future.

Your second point leads to two options:
i) providing face integration functions, but checking inside for evaluate_face_integrals() or has_interior_face_integrals similar as we do for checking the boundary type (Dirichlet, DirichletCached, etc.)
ii) or checking outside and skipping the face integration function argument to matrix_free->loop(). Probably the first would look better and move some if statements inside ... but I think that is something that would need some measuring, or maybe you have done that already in the past, since you accepted the current state of the ifs being inside/close to some kernels.
Here, I cannot guess how heavy the impact actually is, but I suggest renaming as described above and would be willing to do some tests with @kronbichler regarding the second point!

[kronbichler]

This is a somewhat delicate topic. The problem is that inner faces create additional dependencies, both for the ghost value update and the scheduling of operations before/after the loop over cells. In dealii/dealii#15381 I fixed one of the issues, namely the fact that inner face were always generated, also for cases when one did not even request them.

@kronbichler is it a measureable performance penalty when using matrix_free->loop() but leaving e.g. the function with the face integrals empty?

It is definitely measurable, especially for the more advanced algorithms. For the simple case of cell_loop versus loop, the former can be more aggressive in scheduling the pre/post loops around the data access, because the cells have a much smaller domain of dependence. The former also significantly reduces the data accessed by ghost value operations. So it certainly does pay off to put an if into the single base operator instance. I have not entirely made up my mind yet how a good scheduling for the combination of cell integrals + boundary integrals should look like, but I think we should definitely try to address it. One could reason whether it is possible to solve this entirely in deal.II without user-facing differences also for the non-conforming/conforming case, so I would be happy to have those discussions.

In that sense, checking inside evaluate_face_integrals() is definitely too late for deal.II to handle the case, because then it already switches to the conservative algorithms that are potentially considerably slower. So the case (ii) mentioned by @richardschu is the better option.