Implement the Fast Diagonalization Method for Q_k, DQ_k and RTCF_k (FDMPC)

[pbrubeck]

FDMPC is a preconditioner for tensor-product elements that changes the shape functions so that the H^1 Riesz map is diagonalized in the interior of a Cartesian cell, and assembles a global sparse matrix on which other preconditioners, such as ASMStarPC, can be applied.

Other relevant fixes and features coming with this PR:

• Added ASMExtrudedStarPC: works for uniform extrusions, variable layers are not supported.
• pmg.py caused import numpy as np on every script which had from firedrake import * along with many other imports. These imports are dropped by specifying __all__ in pmg.py.
• Enhanced p-coarsening of finite elements: by default for MixedElements and TensorElements, we extract the highest degree and reconstruct the sub-elements while preserving the differences of their degrees relative to the highest degree. This enables us to coarsen spaces like NCF(N) x DQ(N-1).
• BLAS transfer kernels extended to tensor product H(div) and H(curl) elements.
• p-coarsening of the nullspace: here the right thing to do is to restrict (as opposed to inject) the nullspace basis. If we want a solution orthogonal to the nullspace in the L2 sense, the nullspace basis vector must be provided as a co-function which might differ from the nullspace of the discrete operator. This will not occur if we provide constant=True as the basis.
• StandaloneInterpolationMatrix and MixedInterpolationMatrix were defined as restriction matrices before this PR.

[pbrubeck]

Removing the numpy import will potentially affect a lot of users, but this should be very easy to fix. The application libraries can be imported, but numpy might potentially be missing in some of their tests.

[wence-]

Removing the numpy import will potentially affect a lot of users, but this should be very easy to fix. The application libraries can be imported, but numpy might potentially be missing in some of their tests.

I am confused by this:

from firedrake import *

numpy => NameError

At least on my firedrake.

[wence-]

Removing the numpy import will potentially affect a lot of users, but this should be very easy to fix. The application libraries can be imported, but numpy might potentially be missing in some of their tests.

I am confused by this:

from firedrake import *

numpy => NameError

At least on my firedrake.

Ah, np, dropping that is not a problem.

[wence-]

I had mostly minor comments.

I do not claim that I fully understand all the intricacies in this code. I think the core idea is the following:

1. There exists a basis for Q_k (etc...) in which certain bilinear forms are actually sparse (as sparse as the Q_1 thing approximately, so we can assemble them as sparse matrices)
2. This preconditioner works by applying the change of basis (which can be done dimension-by-dimension in a O(p^{d+1}) way from standard GLL to this special basis, call it T_sp
3. In the special basis you assemble the operator taking care to maintain the sparsity, this can now be sent to PETSc (call this operator A_sp)

So the action of the preconditioner is something like:

y <- T_sp^T A_sp^-1 T_sp x

You have to do a lot of work by hand here to construct A_sp correctly (mostly) including making certain approximations on coefficients in the form. I think (though do not know) that if you implemented the fancy 1D basis in FIAT, it should be possible to deliver T_sp as the interpolation between two finat elements (and this should now have good complexity).

The main work in A_sp is that you're having to do all of the assembly by hand, because we can't deliver you an element kernel that maintains sparsity, and even if we could do that, can't deliver the insertion in a sparsity-preserving way.

It would be interesting (not for this PR) to think about what would be required to do this.

[pbrubeck]

You have to do a lot of work by hand here to construct A_sp correctly (mostly) including making certain approximations on coefficients in the form. I think (though do not know) that if you implemented the fancy 1D basis in FIAT, it should be possible to deliver T_sp as the interpolation between two finat elements (and this should now have good complexity).

Yes, I would be very keen to implement this idea. Furthermore we need to figure out a way to take the transpose, which would be very useful for p-MG.

[dham]

Very minor issues.

Addressed