When discretising a PDE using most variants of the spectral element method, the resulting problem is usually expressed as a matrix equation. In traditional linear finite element codes, the matrix is usually represented as a large sparse global matrix, which represents the action of a particular operator such as the Laplacian matrix across the whole domain.
However, when we consider spectral element methods, in which the polynomial order representing the expansion can be far higher, this method becomes far less optimal. We can instead consider the action of an operator locally on each element, and then perform an assembly operation. This is mathematically equivalent to the global matrix approach and gives exactly the same answer, but at high polynomial orders it is far more efficient on modern CPU architectures.
Furthermore, this local approach can be represented in one of two ways: either as a dense matrix for each element, which is typically more efficient at intermediate polynomial orders, or in the hp element case as a tensor product of smaller dense matrices via an approach deemed sum-factorisation, which is used at very high polynomial orders. Figure ?? gives an overview of these three different operator strategies.
A goal of Nektar++ is to support not only high order expansions, but all orders from low (where element size h is the dominant factor) to high (where p dominates); a procedure we have dubbed “from h to p efficiently”.
An obvious question is: “which strategy should I select?” The most important factors in this decision are:
Generally you can use results from three publications [33, 6, 5] which outline results for two- and three-dimensional elements.
In general, the best approach is to perform some preliminary timings by changing the appropriate variables in the session file, which is outlined below. As a very rough guide, for 1 ≤ p ≤ 2 you should use the global approach; for 3 ≤ p ≤ 7 use the local approach; and for p ≥ 8 use sum-factorisation. However, these guidelines will vary due to the parameters noted above. In future releases of Nektar++ we hope to tune these variables automatically to make this decision easier to make.
Operator evaluation strategies can be configured in the GLOBALOPTIMISATIONPARAMETERS
tag,
which lies inside the root NEKTAR
tag:
Operator evaluation is supported for four operators: backward transform, inner product, Helmholtz and mass operators. It is possible to specify the following optimisation flags for different operators:
DO_GLOBAL_MAT_OP
: If VALUE
is 1
, the globally assembled system matrix will
be used to evaluate the operator. If VALUE
is 0
, the operator will be evaluated
elementally.
DO_BLOCK_MAT_OP
: If VALUE
is 1
, the elemental evaluation will be done using
the elemental/local matrices (which are all concatenated in a block matrix,
hence the name). If VALUE
is 0
, the elemental evaluation will be done using the
sum-factorisation technique.
Each element type (triangle, quadrilateral, etc) has its own VALUE
, since break-even
points for sum-factorisation and the local matrix approach will differ depending
on element type. Note that due to a small shortcoming in the code, all element
types must be defined; so three-dimensional elements must still be defined even if
the simulation is two-dimensional.
Note that global takes precendence over block, so if VALUE
is set to 1
for both then the
operator will be global.
For very complex operators – in particular HelmholtzMatrixOp
– always set DO_BLOCK_MAT_OP
to
1
as sum-factorisation for these operator types can be costly.