Improve compute penalty parameter

[kronbichler]

In some FSI cases, @nfehn and I observed very high cost of initializing FEFaceValues (for high-order mappings). This PR avoids this by computing the volume and surface area by the much faster FEEvaluation and FEFaceEvaluation concepts. The difficulty is that we need to exchange data between different MPI processes, for which this PR suggests to use the global active cell index (avaiable since the 9.3 release of deal.II) and a distributed vector. This should avoid the bottleneck entirely and make the function as fast as a cell loop.

[nfehn]

I only have (several) minor comments. Maybe also @peterrum could have a brief look at this PR once you have addressed the requested changes.

[peterrum]

I'll do. @kronbichler could you ping me once you are ready?

[kronbichler]

In 3449da6 I added better comments and re-organized the code: I now first compute the surface areas and then go through the cells, which allows me to directly store the result back without keeping the intermediate volume in the vector.

[nfehn]

@peterrum Here is the ping :)

[kronbichler]

I can make those changes, but before we spend more time on this side topic, let me raise a question: What would you think if we introduce a special parameter for quad_index in FEEvaluation (either the "default" value or a special one we can provide), https://github.com/dealii/dealii/blob/be72dbb5a5f97aee1c8cf330b66664c848176fba/include/deal.II/matrix_free/fe_evaluation.h#L2049, that picks up the quadrature formula which is closest (from above) to degree + 1? This might be numbers::invalid_unsigned_int (similar to the hp parameters) or a more expressive name, e.g. FEEvaluation::detect_quadrature_formula_from_degree. I am not trying to hide from this discussion here, but I would like to think about the best solution for the long term and if we can make usability improvements, given this experience.

[peterrum]

I can make those changes, but before we spend more time on this side topic, let me raise a question: What would you think if we introduce a special parameter for quad_index in FEEvaluation (either the "default" value or a special one we can provide), https://github.com/dealii/dealii/blob/be72dbb5a5f97aee1c8cf330b66664c848176fba/include/deal.II/matrix_free/fe_evaluation.h#L2049, that picks up the quadrature formula which is closest (from above) to degree + 1? This might be numbers::invalid_unsigned_int (similar to the hp parameters) or a more expressive name, e.g. FEEvaluation::detect_quadrature_formula_from_degree. I am not trying to hide from this discussion here, but I would like to think about the best solution for the long term and if we can make usability improvements, given this experience.

I was thinking about the same yesterday. We also had something similar at one time in the code (I don't recall in which context): we looped over all quadrature rules and checked the number of quadrature points. The problem I saw were:

• how does the generalization look like: simplex, wedge, pyramid, ...
• let's assume we have two sets of quadrature rules for degree + 1: one with Gauss–Legendre quadrature and the other with Gauss–Lobatto quadrature. Which one do we want to take? Furthermore, the result will depend on the order how the quadrature rules have been registered. Probably, choosing any of the two quadrature rules will work but might lead to confusions in the future...

[nfehn]

I am not sure whether more automation provides better usability in this context. Personally, I don't like initializing variables with some special numbers or invalid numbers in order to express that a special path should be taken. The code becomes difficult to understand and, most often, it requires a dictionary for translation.

[kronbichler]

In the constructor of FEEvaluation we can only afford a reasonably cheap operation. A solution actually looking into the quadrature formulas and determining whether it is the most accurate would not work. So this brings us back to the other discussion we had in terms of the overall evaluator design, allowing a user to collect all information to attach to a particular evaluation context and not just plain integers.

[nfehn]

I suggest to complete this PR first and get is merged soon. Only a few additional lines of code should be necessary as described in detail in the code review.

[kronbichler]

I rebased this PR and added the quad_index_standard to the convection-diffusion and incompressible Navier-Stokes modules. I am not very happy about the fact that we have to pull this parameter through so many interfaces, given that it is only for computing a parameter based on a rough estimate. More precisely, I am comparing the state after the last (fourth) commit to the one after the second commit (where we simply pick the zero-th quadrature index, and take the possibility of a small quadrature error into account, modulo some interface cleanup that I addressed in the third commit during rebasing). At least it is clearly better than after the third commit, so we might go on from here.

[nfehn]

I agree with your thoughts on pulling the quad index through so many interfaces, and that we would like to avoid this. I think the solution is adding it to the respective KernelData. This way, the combined operators for convection-diffusion and Navier-Stokes should have all the information they need, given the modular design that we have chosen for these combined operators (through composition used in the Operator and the OperatorData classes).

[nfehn]

@kronbichler Is it ok for you that we try to merge #296 first, in order to make progress with simplex?