Low-Order Refined Preconditioning for Spectral/hp Element Method for Complex, 3D Geometries

1. Introduction.

• The preconditioner constructs a refined low-order P=1 finite element discretisation that is spectrally equivalent to the higher-order system.The paper implements this approach within Nektar++ for incompressible Navier–Stokes problems.Source: §1 · p. 1
• LOR uses the sparsity of low-order operators to build a sparser preconditioner, with constant memory per degree of freedom and reduced memory use compared with higher-order discretisations.The cited prior work reports bounded iterative condition-number behaviour under specific conditions and minimal sensitivity to high-aspect-ratio elements.Source: §1 · p. 2
• Earlier LOR developments focused especially on rectangular and hexahedral meshes, including work on GPU acceleration and scalability [31].The cited work also covers H(curl) and H(div) finite element problems using N´ed´elec and Raviart–Thomas elements.Source: §1 · p. 2
• At high polynomial order, GLL point spacing becomes strongly non-uniform, creating anisotropic elements that challenge standard multigrid and domain-decomposition methods.The spacing scales as 1/P^2 near element boundaries and 1/P in the interior, and solver efficiency can deteriorate on spectrally refined meshes.Source: §1 · p. 2
• The implementation supports mixed two- and three-dimensional meshes, including simplicial elements, and extends LOR to prismatic and tetrahedral meshes and both Lagrange and hierarchical modal bases.The framework is designed for large-scale computations and supports problem sizes exceeding 10^7 degrees of freedom.Source: §1 · p. 3

2.1. Incompressible Navier-Stokes equations.

• The velocity-correction scheme advances incompressible flow through explicit convection, an implicit pressure Poisson solve, and implicit Helmholtz solves for velocity components.The pressure is an auxiliary scalar used to enforce the divergence-free velocity condition.Source: §2.1 · p. 3

2.2. The spectral/hp element method.

• The spectral/hp method combines finite-element geometric flexibility with spectral convergence through polynomial-order basis functions on element-partitioned domains.The work uses a hierarchical modified C0 basis and Gaussian or GLL point distributions for integration.Source: §2.2 · pp. 3–4

2.3. Solving the linear system.

• The pressure Poisson system is elliptic and ill-conditioned on large, complex anisotropic meshes, motivating preconditioning to reduce the modified system’s condition number.The velocity-correction scheme solves an implicit pressure Poisson problem to enforce incompressibility.Source: §2.1 · p. 3, §2.3 · p. 4

2.4. LOR preconditioner.

• The LOR construction replaces dense within-element high-order coupling with connections between immediate neighbours, making its operators cheaper to evaluate.The refined grid places P=1 elements at higher-order quadrature points and preserves their GLL-based point distribution.Source: §2.4 · p. 5, Figure 1 · p. 5

2.4.1. Interpolation to LOR space.

• The high-order modal space is mapped to the LOR space through a Vandermonde basis transformation followed by collocation at shared GLL nodes.At the collocated nodes, the high-order nodal and piecewise-linear LOR coefficient vectors coincide.Source: §2.4.1 · pp. 5–7

2.4.2. Operator Equivalence and Preconditioning.

• The modal and nodal preconditioned operators have identical spectra, so Krylov convergence is independent of the high-order basis.Any practical difference is attributable to the conditioning of the Vandermonde matrix; preconditioner construction is a one-time Build stage whose costs are amortised over iterative solves.Source: §2.4.2 · p. 8, §3 · p. 8

3. Algorithm.

• Each preconditioner application maps the high-order residual into LOR space, solves the assembled LOR system, then projects and gathers the correction back to the high-order space.Multiplicity corrections account for shared degrees of freedom, particularly the repeated mappings in simplicial meshes.Source: §3 · pp. 8–10
• The one-time build stage constructs and caches the refined mesh, mappings, LOR system, and preconditioner before repeated application.For a modal high-order basis, construction also requires the generalised Vandermonde matrix.Source: §3 · p. 10

3.1.1. LOR space for 2D element types.

• For simplices, the LOR construction uses electrostatic interior points and boundary points constrained to neighbouring faces or edges; prisms combine triangular and one-dimensional GLL subdivisions.The triangular nodal expansion is combined with a one-dimensional Lagrange expansion for prisms.Source: §3.1 · p. 11, Figure 6 · p. 12, §3.1.2 · p. 13

3.1. Choice of LOR space.

• Quadrilateral LOR refinement offers Q1 quadrilateral or P1 triangular subelements on the GLL points; triangles are subdivided only into smaller P1 elements.The P1 quadrilateral option produces 2n^2 triangular subelements.Source: Figure 4 · p. 11, §3.1.1 · pp. 12–13

3.1.2. LOR space for 3D element types.

• The inner solver can use a direct method for small cases or a single MPI rank, while larger cases can use PETSc or native Nektar++ solvers and preconditioners.Available PETSc options include BoomerAMG and GAMG.Source: §3.2 · p. 13

3.2. Inner-iteration strategy.

• The inner-iteration strategy affects LOR-preconditioner performance, and the pressure Poisson system is the study's test problem.The pressure Poisson system is identified as the primary performance bottleneck in the incompressible Navier–Stokes solver.Source: §3.2 · p. 13, §4 · p. 13

4.1.1. 2D test cases.

• For 2D meshes, triangular condition numbers increase exponentially and exceed quadrilateral results at P ≥5, while quadrilateral maximum eigenvalues remain controlled and triangular maxima increase.The minimum eigenvalue stays constant across the tested cases; stronger high-order LOR-space anisotropy may make high-frequency components harder to condition for triangles.Source: §4.1.1 · pp. 14–15, Figure 8 · p. 15

4.1.2. 3D test cases.

• In 3D eigenvalue distributions, the prism minimum stays approximately constant while the tetrahedron minimum falls and its maximum rises sharply as P increases.For both element types, eigenvalue clustering density remains similar while the maximum eigenvalue increases.Source: Figure 10 · p. 16, §4.1.2 · p. 16
• On prism and tetrahedron cube meshes, LOR condition numbers are more than an order of magnitude below diagonal preconditioning and remain controlled through P ≤6.The gap is particularly notable for P ∈[3, 6].Source: §4.1.2 · p. 16

4.2.1. 2D NACA aerofoil.

• At P = 4, LOR requires the fewest iterations for the NACA0012 case, while mixed-mesh iteration counts remain controlled despite losing bounded growth with polynomial order.The loss of bounded behavior is attributed to triangles lacking an optimum LOR-space node set.Source: Table 2 · p. 17, §4.2.1 · pp. 17–18

4.2.2. 3D Hemisphere.

• For P ≤6, LOR has the best iterative performance across both hemisphere meshes and boundary-condition choices.At higher orders, diagonal-preconditioned GMRES performs better than LOR only for the Dirichlet case.Source: Figure 13 · p. 18, §4.2.2 · p. 19
• Using one AMG V-cycle per LOR inner iteration gives outer iteration counts not far from solving the inner system to a 10^-4 residual drop, though counts vary across runs.The tighter residual uses multiple V-cycles and yields fewer outer iterations.Source: §4.2.2 · p. 19

5.1. Strong scaling performance.

• Outer CG iteration counts are essentially independent of core count at P = 2 and P = 3, carrying the established conditioning behaviour to industrial scale.The strong-scaling sweep starts at 512 cores for P = 3 because 256-core runs exhausted available memory.Source: §5.1 · pp. 20–21
• The one-time AMG setup grows from 12s to 197s for P = 2 and from 43s to 233s for P = 3 as core count increases, but is reused across production timesteps.The authors exclude setup from time-to-solution and efficiency metrics because the pressure operator is constant-coefficient in production.Source: §5.1 · p. 21
• For P = 2, KSPSolve time falls from 14.3s at 256 cores to 4.7s at 2048 cores; for P = 3, it falls from 40.1s at 512 cores to 16.3s at 4096 cores.Efficiency at 2048 cores is 38.2% for P = 2 and 59.9% for P = 3; solve-time metrics exclude setup cost.Source: Figure 15 · p. 20, §5.1 · p. 21

5.2. Incompressible Navier-Stokes solver performance.

• 11.2× lower pressure-solve wall-clock time with LOR than Baseline-full accompanies a reduction in mean outer iterations from 559.6 to 5.5.Baseline-full spends 74% of execution time in the pressure solve, compared with 17% with LOR.Source: §5.2 · p. 21, Figure 16 · p. 22
• 235.3 to 5.5 mean outer iterations with LOR versus Baseline yields only a 1.2× pressure-solve speedup: 2727 s to 2289 s.Each LOR iteration costs 0.41 s versus 0.012 s for Baseline, offsetting most of the iteration-count advantage.Source: §5.2 · p. 22
• Against Baseline, LOR reduces mean pressure iterations by a factor of 42, while static condensation and LOR/AMG address different computational costs.The three configurations average 69.7 velocity iterations, indicating the pressure-preconditioner comparison does not alter velocity-solve convergence.Source: §5.2 · pp. 22–23
• Initial-guess extrapolation across timesteps leaves only the first few of 1000 solves near the 54-iteration cold-start cost, against a production mean of 5.5.The one-time AMG build costs 130s in production versus 136s in the isolated P = 3 benchmark; the effective LOR cost is 2.3s per step.Source: §5.2 · p. 23
• The formulation relates modal and nodal preconditioned operators by a similarity transformation, giving identical spectra and algebraically equivalent preconditioned CG iterations in exact arithmetic.The result is independent of element type, including types without an available spectral-equivalence bound.Source: §6 · p. 23

Conclusions.

• On quadrilaterals, iteration counts show the asymptotically optimal behaviour predicted by tensor-product theory; on triangles, iteration growth remains controlled despite increasing condition number.In 3D, counts remain controlled up to P = 6 on tetrahedral, prismatic, and mixed-element meshes, with fewer iterations than the default Nektar++ configuration.Source: §6 · p. 23
• On the IFW-W mesh of 2.87 × 10^6 mixed prismatic and tetrahedral elements, LOR reduces mean pressure CG iterations at P = 3 from 235.3 to 5.5.Across 1000 production timesteps, pressure-solve time decreases by 16.1% relative to the default preconditioner.Source: §6 · pp. 23–24
• Outer iteration counts remain essentially independent of MPI rank count in strong-scaling tests, while solve-time gains diminish at the highest core counts as communication costs increase.AMG hierarchy construction becomes more expensive with rank count, but is incurred once and amortised over production timesteps.Source: §6 · p. 24
• On simplices, increasing P clusters interpolation points near vertices and edges, producing anisotropic sub-elements associated with increasing condition number.Alternative or optimisation-based node distributions are proposed for further study at higher orders; the current implementation requires an assembled low-order operator.Source: §6 · p. 24

Appendix A. Linear solvers in Nektar++.

• All numerical experiments using the LOR preconditioner were performed with the master branch of Nektar++.Source: Appendix B · p. 24

Appendix B. Software.

• The LOR implementation is provided in Nektar++ merge request !1888.Source: Appendix B · p. 25
• HX1 used PETSc v3.22.0 with HYPRE v2.31.0, while ARCHER2 used PETSc v3.19.3 with HYPRE v2.28.0.Source: Appendix B · p. 25
• The software environment summary covers HX1 and ARCHER2, the HPC platforms used for the studies.Source: Appendix C · p. 25, Table 4 · p. 25

• HX1 nodes have 64 cores, 512 GB of RAM, and a 200 Gbit/s InfiniBand interconnect.Each node has two 32-core Intel Xeon Platinum 8358 processors running at 2.60 GHz.Source: Appendix C · p. 25
• HX1 provides 64 cores and 512 GB per node, while ARCHER2 provides 128 cores and 256 GB per node.The clusters also use different interconnects: NVIDIA ConnectX-6 HDR200 on HX1 and HPE Cray Slingshot on ARCHER2.Source: Appendix C · p. 25
• ARCHER2 nodes use two 64-core AMD EPYC 7742 processors, 256 GB of RAM, and a 2 × 100 Gbps bidirectional Slingshot interconnect.All processor cores are fully utilised during the runs.Source: Appendix C · p. 25

• AMG uses one V-cycle, HMIS coarsening, a 0.70 strong threshold for 3D cases, and a coarse-grid size of 10.The coarse-grid size is set to control the number of AMG levels.Source: Appendix D · p. 25
• Interpolation uses Extended+i with Pmax = 4 and truncation of 0.3.Source: Appendix D · p. 25
• Smoothing uses L1Scaled SOR/Jacobi with two sweeps in both the upward and downward directions.Source: Appendix D · p. 25

• The wind-tunnel domain uses a uniform Dirichlet inlet, slip ceiling and side walls, a stable high-order zero-Neumann velocity outlet, and ambient pressure.The simulation setup follows the methodology detailed in.Source: Appendix E · p. 25
• The simulations use semi-implicit time stepping and a second-order time-integration scheme with the maximum constant stable time-step.Source: Appendix E · pp. 25–26
• The Taylor–Hood approximation uses velocity order P+1 and pressure order P; simulations run at P=4 and P=3.Initialisation uses a velocity field from a restart solution and sets pressure to 0.Source: Appendix E · p. 26

• For a degree-P parent expansion, P+1 edge points define P refined segments, while shared-edge and face identifiers and orientations are preserved.The construction creates required interior entities and updates coefficient maps for transferring data between high-order and LOR representations.Source: Appendix F · p. 26
• The two-dimensional construction subdivides triangular faces using prescribed nodal points and refines quadrilaterals into either two triangles or one quadrilateral per refined cell.It updates local-to-global connectivity and the two-dimensional coefficient map.Source: Appendix F · p. 27
• Tetrahedral subdivision evaluates physical coordinates at simplex nodal points on edges, faces, and interiors, reuses stored face connectivity, and assembles linear tetrahedra.The algorithm updates the global mesh and three-dimensional coefficient map.Source: Appendix F · p. 28
• Prismatic subdivision checks triangular-face alignment, refines edges and faces, and connects corresponding triangular subdivisions into linear prisms with inherited orientations.The algorithm records new offsets and updates the global mesh and geometric connectivity.Source: Appendix F · p. 29