Block jacobi preconditioner

[bergbauer]

This PR makes the BlockJacobiPreconditioner more consistent with the JacobiPreconditioner and the AMG preconditioners. The block matrices should not be owned by the OperatorBase but by the Preconditioner.

The second commit makes separate loops over cells and faces possible in parallel. We should discuss if we want to enable this option as it uses some knowledge about internal datastructures in MatrixFree.

A marked this as a draft because I want some feedback before I put more work into this.

[nfehn]

Principally, I like the idea that the preconditioner takes care of the data, just as the Jacobi preconditioner takes care of the diagonal. So see whether this design has downsides, I need to go through the changes in detail. I will add questions as comments to the code changes.

[nfehn]

Here, it looks like the BlockJacobiPreconditioner always uses a matrix-based implementation technique. In case a matrix-free application of the inverse block "matrices" is chosen, block_matrix is probably not used. We should make sure how to make this intuitively understandable when looking at the implementation of BlockJacobiPreconditioner.

[nfehn]

Shouldn't this function, and the other one below for the inverse block matrices, be implemented as a free function? So we need to discuss: Why is OperatorBase executing this function? With the changes of this PR, isn't it unsatisfactory that this is still done by OperatorBase?

[nfehn]

here, a similar comment applies as for the BlockJacobiPreconditioner.

[nfehn]

I think the present PR moves the matrix-based and the matrix-free implementation of Block-Jacobi operations further away from each other than they are currently. From the matrix-based perspective, storing the block matrices outside the operator makes definitely sense. From the matrix-free perspective, the operations are very close to what OperatorBase can do, due to the tight coupling to MatrixFree. We will probably not be able to remove block-Jacobi functionality from OperatorBase (meaning computing the relevant integrals in a matrix-free way). The question then is: do we consider block-Jacobi operations as a black-box linear algebra functionality, or as a simplified version of the discrete PDE operator (just the block-Jacobi part of the whole operator)?

A difference to the Jacobi preconditioner is that we always store the diagonal, even if we use a matrix-free implementation framework. This is different for block-Jacobi.

Currently, I ask myself whether it makes sense to unify the matrix-based and matrix-free Block-Jacobi operations in a common data structure, say BlockJacobiOperator. This data structure would contain the block-matrices, but also all the Elementwise data structures needed for the matrix-free implementation of block-Jacobi. The question is whether we can reasonably pull out the matrix-free case out of OperatorBase, or whether the coupling to OperatorBase/MatrixFree is too tight? In the present PR, we should definitely elaborate on this and give this a try. Ideally, OperatorBase only provides functionality to perform different types of evaluations of operators and integrals (including block-Jacobi operations), but the data we operate on is living outside OperatorBase, so no block-Jacobi matrices and also no elementwise operator/preconditioner/solver data structures in OperatorBase. According to such a design, BlockJacobiOperator would receive an instance of OperatorBase as constructor argument, in order to be able to perform the relevant integrals. Multigrid or the block-Jacobi preconditioner would only work with / see a BlockJacobiOperator (which can be matrix-based or matrix-free; multigrid or the preconditioner need not know this).