Need clarification on a piece of book excerpt about spectral element method!

Question

I am reading "Using MPI (3rd edition)" from William Gropp, where in chapter 4 application section 4.13, it introduces an MPI application Nek5000/NekCEM which is based on spectral element method (SEM) coupled with high-order timestepping to solve the governing IVM/BVM problems.

There is one statement about SEM that I do not quite understand (Yes, it's not relevent to MPI):

"A central idea in the SEM is never to form local elemental matrices, which would in general have $(N+1)^6$ nonzeros per element, but to instead use preconditioned iterative solvers that requrie only the action of operators applied to functions."

I have some limited backgroud in FEM and SEM, and I know that SEM elemental matrices are usually much larger than FEM due to high polynomial orders, but I still can't get what the author is trying to conveying with this statement.

Could someone please explain it with more concrete descriptions? Any help would be greatly appreciated!

Answer 1

In SEM there is usually an explicit formula for the action of the local matrices on local vectors. The local matrices can be thought of as numerically precomputed operators that evaluate: interpolation, differentiation, and integration. People precompute them in tabulated form as matrices because that is the only way to produce a sparse matrix, which is the format direct solvers like LU,Cholesky,... expect. For SEM though these local matrices require a lot of memory, and iterative solvers do not require a sparse matrix input, so we can avoid precomputing these operators.

One way to avoid precomputing these is to compute interpolants, integrals, and derivatives by way of three-term-recurrences from underlying orthogonal polynomial basis. In SEM most basis have such a recurrence relation, or are lightly modified versions of a basis which does.

In many cases also you only need to store "reference" matrices by precomputing these operators on a reference element, and then compute the desired actions of these by geometric transformations alone. You can store geometric transformation factors for each element of the mesh, or even recompute that information at need to save even more memory.

This only does the local action though. It's still up to you to get inter-element details right. You will need to store some kind of connectivity graph(s) for your mesh and use that information to piece together all the local matrix-vector products into a full action of the operator.

That full operator evaluation can be fed into a Krylov solver as a black-box function, and that is enough information for it to try and solve your resulting system.