Fact-checked by Grok 2 weeks ago

Spectral element method

The spectral element method (SEM) is a high-order numerical technique for solving partial differential equations that integrates the geometric flexibility of the finite element method with the exponential accuracy of spectral methods, employing high-degree piecewise polynomial basis functions defined over a mesh of elements.^[1] This approach uses nodal basis functions, such as Lagrange polynomials of degree 4 to 10, evaluated at Gauss-Lobatto-Legendre points to achieve precise approximations with fewer degrees of freedom compared to low-order methods.^[2] Originally developed in the field of computational fluid dynamics in 1984 by Anthony T. Patera, the SEM was later adapted for seismology in the 1990s to model wave propagation in complex three-dimensional Earth structures.^[2] It employs a weak variational formulation of the governing equations, often resulting in a diagonal mass matrix through exact quadrature, which enhances computational efficiency and stability.^[2] Variants like the spectral/hp element method further allow adaptive refinement in both mesh size (h) and polynomial order (p), enabling exponential convergence for smooth solutions under sufficient regularity.^[3] Key advantages of the SEM include superior accuracy for problems involving smooth solutions, reduced numerical dispersion and dissipation, and the ability to handle irregular geometries, material heterogeneities, and interfaces such as fluid-solid boundaries without iterative coupling.^[4] These features make it particularly effective for large-scale simulations on parallel architectures, as demonstrated in global seismic modeling with billions of grid points.^[2] Applications span seismology for synthetic seismogram generation, computational fluid dynamics for transitional flows, structural dynamics in soil-structure interactions, and electromagnetic problems like magnetotelluric modeling.^[4]

Overview

Definition and Principles

The spectral element method (SEM) was first proposed by A. T. Patera in 1984 as a hybrid numerical technique extending ideas from both spectral methods and finite element methods, initially applied to fluid dynamics problems such as laminar flow in channel expansions.^[5] It serves as a high-order variant of the finite element method, utilizing piecewise polynomials of degree N within each element to approximate solutions, thereby merging the geometric flexibility of finite elements for handling complex domains with the superior accuracy of spectral techniques for smooth functions.^[6] At its core, SEM relies on principles of high-order approximation, where the use of elevated polynomial degrees enables spectral-like accuracy—characterized by exponential convergence rates—locally within elements for sufficiently smooth solutions.^[7] This is complemented by a tensor-product structure in the basis representation, which facilitates efficient tensor-based operations and reduces computational costs, especially in structured or semi-structured domains like those encountered in computational fluid dynamics or wave propagation.^[8] In practice, SEM approximates solutions to partial differential equations (PDEs), including the Navier-Stokes equations governing incompressible flows, by employing variational formulations that project the continuous problem onto a discrete polynomial subspace, ensuring conformity and stability through weighted residuals.^[5] This approach allows for robust handling of both elliptic and time-dependent problems while maintaining high fidelity in regions of interest.^[6]

Advantages and Limitations

The spectral element method (SEM) exhibits exponential convergence rates for smooth solutions, achieving error reductions on the order of O(\exp(-\gamma N)), where N is the polynomial degree and \gamma > 0 is a constant depending on the solution regularity.^[9] This high-order accuracy allows SEM to attain a target precision with significantly fewer degrees of freedom compared to low-order finite element methods (FEM), which typically require O(h^{-d}) degrees of freedom in d dimensions for mesh size h.^[1] In SEM, the degrees of freedom scale as O(N^{d}) per element due to the tensor-product structure of high-order polynomials, enabling efficient solutions for problems where a modest number of elements suffices.^[10] Additionally, the use of Gauss-Lobatto-Legendre quadrature in SEM yields exactly diagonal mass matrices, facilitating fast explicit time-stepping schemes by avoiding matrix inversions or factorizations.^[8] SEM is particularly well-suited for simulations involving smooth solutions and geometries of moderate complexity, such as channel flows in fluid dynamics or layered media in wave propagation, where its spectral accuracy combines effectively with finite element flexibility.^[9] In these contexts, the method delivers high-fidelity results with reduced computational overhead relative to traditional approaches.^[11] Despite these strengths, SEM performs poorly on highly irregular geometries, often necessitating a proliferation of elements to maintain accuracy, which erodes its efficiency advantages.^[12] The higher computational cost per element arises from the tensor-product operations inherent to high-order basis evaluations, leading to denser local matrices and increased expense for large polynomial degrees.^[9] Furthermore, SEM's performance is sensitive to the choice of polynomial degree N, as suboptimal selection can amplify costs without commensurate accuracy gains, particularly in nonsmooth problems.^[1]

Mathematical Formulation

Basis Functions and Polynomial Spaces

In the spectral element method (SEM), the approximation space within each element is constructed using high-degree nodal Lagrange polynomials interpolated at Gauss-Lobatto-Legendre (GLL) points. These points, which include the endpoints of the reference interval, serve as both interpolation nodes and quadrature points, enabling exact numerical integration of the mass matrix for polynomials up to degree 2N when using (N+1) points. This choice ensures a diagonal mass matrix, which significantly simplifies the solution process by avoiding the need to invert dense matrices explicitly.^[5]^[2] The one-dimensional Lagrange basis function centered at the k-th GLL point \xi_k is defined as

l_k(\xi) = \prod_{j=0, j \neq k}^{N} \frac{\xi - \xi_j}{\xi_k - \xi_j},

where \{\xi_j\}_{j=0}^{N} are the roots of the equation (1 - \xi^2) P_N'(\xi) = 0, with P_N(\xi) denoting the N-th Legendre polynomial. These basis functions form a nodal representation, allowing direct evaluation of the solution at the GLL points without additional interpolation. The underlying Legendre polynomials provide orthogonality properties that facilitate efficient computation of derivatives through spectral differentiation matrices, which approximate spatial derivatives exactly for polynomials of degree up to N.^[2] In two or three dimensions, the basis functions exploit a tensor-product structure, where the multivariate basis is the product of univariate Lagrange polynomials along each coordinate direction within hexahedral elements. This results in (N+1)^d nodes per element, with d being the spatial dimension, enabling high-order accuracy while maintaining geometric flexibility. Unlike modal bases, which expand solutions in terms of orthogonal polynomials (e.g., Legendre or Chebyshev modes) and often require iterative solvers for the mass matrix, the nodal Lagrange basis in SEM supports straightforward pointwise evaluation and is particularly advantageous for problems with non-uniform meshes or complex geometries.^[5]

Galerkin Discretization and Weak Form

The spectral element method (SEM) applies the Galerkin projection to obtain a discrete variational formulation of partial differential equations (PDEs). Consider the model Poisson equation -\Delta u = f in a domain \Omega \subset \mathbb{R}^d with suitable boundary conditions. The weak formulation requires finding u \in H^1(\Omega) such that

\int_{\Omega} \nabla u \cdot \nabla v \, dx = \int_{\Omega} f v \, dx

for all test functions v \in H^1(\Omega).^[13] This variational form shifts derivatives from the solution to the test functions, enabling the use of finite-dimensional subspaces for approximation. In SEM, the solution is approximated in a finite element space V_N^h consisting of continuous, piecewise polynomials of degree at most N on a partition of \Omega into E non-overlapping spectral elements \{\Omega_e\}_{e=1}^E. The SEM solution u_N \in V_N^h satisfies the discrete weak form

\int_{\Omega} \nabla u_N \cdot \nabla v \, dx = \int_{\Omega} f v \, dx \quad \forall v \in V_N^h,

known as Galerkin orthogonality, which ensures that the residual is orthogonal to the approximation space. The approximate solution takes the nodal form u_N = \sum_i u_i \phi_i, where \{\phi_i\} are the global basis functions spanning V_N^h (as defined in the section on basis functions and polynomial spaces). The integrals in the weak form are evaluated element-wise to assemble the global system. On each element \Omega_e, the local stiffness matrix entries are A^e_{ij} = \int_{\Omega_e} \nabla \phi_i \cdot \nabla \phi_j \, dx and, for mass terms in more general PDEs, M^e_{ij} = \int_{\Omega_e} \phi_i \phi_j \, dx. These are computed using Gauss-Lobatto-Legendre (GLL) quadrature with (N+1)^d points per element in d dimensions, which integrates products of polynomials of degree up to $2N exactly due to the properties of the Lagrange basis at GLL nodes.^[2] The global stiffness matrix is then formed by summing contributions from all elements, respecting C^0 continuity at element interfaces. Boundary conditions are handled distinctly in the weak form. Essential (Dirichlet) conditions are enforced strongly by prescribing nodal values u_i at boundary GLL points, modifying the degrees of freedom accordingly. Natural (Neumann or Robin) conditions appear as boundary integrals in the weak form and are incorporated weakly during assembly without explicit modification of the basis. For time-dependent PDEs, such as the heat equation \partial_t u - \Delta u = f, the SEM extends the spatial Galerkin discretization to a semi-discrete system M \partial_t \mathbf{u}_N + A \mathbf{u}_N = \mathbf{f}, where M is the global mass matrix and A the stiffness matrix. The GLL quadrature yields a diagonal M, facilitating efficient explicit time-stepping schemes like Runge-Kutta methods via simple matrix-vector multiplications, or implicit schemes solved with fast diagonal inversions.^[2] This collocated structure enhances computational efficiency for hyperbolic and parabolic problems.

Numerical Implementation

Domain Discretization and Meshing

In the spectral element method (SEM), the computational domain is discretized into a collection of non-overlapping elements, primarily quadrilateral in two dimensions and hexahedral in three dimensions, to leverage the efficiency of tensor-product structures for high-order polynomial approximations.^[2]^[14] These element types enable separable operations in each coordinate direction, reducing computational cost compared to simplicial elements like triangles or tetrahedra.^[2] Each physical element is mapped from a reference element, typically the cube [-1, 1]^d where d is the spatial dimension, using isoparametric transformations that employ the same basis functions for both geometry and solution representation.^[15]^[16] This mapping preserves the tensor-product form within the reference domain while accommodating curved boundaries through higher-order geometric descriptions.^[15] Meshing strategies in SEM prioritize structured grids for simple geometries, such as rectangular or cubic domains, where elements align with coordinate axes to minimize distortion and facilitate exact tensor-product evaluations.^[17] For domains of moderate complexity, unstructured meshes composed of quadrilateral or hexahedral elements are employed, often with affine mappings to maintain low distortion and ensure numerical stability.^[18] High-distortion elements, which can arise from aggressive curving or poor aspect ratios, are generally avoided to prevent degradation in accuracy, as they introduce errors in the mapping Jacobian that affect integration and conditioning.^[18] In practice, mesh generation tools tailored for SEM, such as those producing body-fitted grids, ensure elements conform to the domain boundaries while keeping the number of elements coarse due to the high polynomial degrees used.^[17] Refinement in SEM typically follows an h-p framework, where a fixed high polynomial degree p (often N \geq 5) is paired with a coarse element size h to achieve exponential convergence, contrasting with the h-refinement dominant in low-order methods.^[19] Adaptive p-refinement allows local increases in polynomial order within elements to resolve features like boundary layers or shocks, while h-refinement adds elements only where necessary, balancing accuracy and cost.^[20] This approach enables efficient resolution on relatively few elements, with examples showing optimal performance using 3x3 meshes for certain problems before varying p.^[21] At element interfaces, C^0 continuity is enforced by sharing Gauss-Lobatto-Legendre (GLL) nodes along boundaries, which serve as both interpolation and quadrature points, ensuring weak coupling without explicit flux computations in continuous Galerkin formulations.^[22]^[23] This nodal matching simplifies the enforcement of inter-element continuity for the solution and its derivatives up to order p-1.^[22] Handling curvilinear elements introduces challenges, as the isoparametric mapping requires computation of the Jacobian determinant for volume integrals and its inverse for gradient terms, elevating setup costs particularly for high p where more quadrature points are involved.^[16]^[24] These computations, performed once per element during preprocessing, can become burdensome for complex geometries but are essential for accurate representation of deformed domains without sacrificing spectral accuracy.^[16]

Assembly, Solving, and Parallelization

In the spectral element method (SEM), the global system is assembled by summing the contributions from local element matrices to form the stiffness matrix \mathbf{K} and mass matrix \mathbf{M}, where the local support of basis functions ensures sparsity in the global matrices. The mass matrix \mathbf{M} is diagonal due to the use of Gauss-Lobatto-Legendre (GLL) quadrature points, which align with the nodal Lagrange basis and exact integration properties for polynomials up to degree $2P-1, where P is the polynomial order. This diagonal structure simplifies operations, particularly in time-stepping schemes, while the stiffness matrix \mathbf{K} remains sparse with bandwidth proportional to the number of elements interfacing each GLL point. To enhance computational efficiency, especially for high-order polynomials, SEM often employs a matrix-free approach that avoids explicit storage of \mathbf{K} and \mathbf{M}, instead computing matrix-vector products on-the-fly using sum-factorization techniques. Sum-factorization exploits the tensor-product structure of element interiors to evaluate integrals via nested 1D operations, reducing the per-element cost from O(P^{3d}) to O(P^{d+1}) in d dimensions, where P approximates the total degrees of freedom N per element. This enables scalable simulations without the memory overhead of dense or sparse matrices, which would scale as O(N^2) globally. For steady-state problems, the resulting linear system \mathbf{K} \mathbf{u} = \mathbf{f} is typically solved using iterative methods such as the conjugate gradient (CG) algorithm, preconditioned with the diagonal of \mathbf{M} or \mathbf{K} to cluster eigenvalues and accelerate convergence. In time-dependent problems, the semi-discrete form \mathbf{M} \frac{d\mathbf{u}}{dt} + \mathbf{K} \mathbf{u} = \mathbf{f} is advanced using explicit Runge-Kutta schemes, where the inverse mass matrix action \mathbf{M}^{-1} is trivial due to its diagonal nature, and operator applications leverage matrix-free sum-factorization for efficiency. These explicit methods support stability up to a CFL-limited time step scaling as O(1/P), balancing accuracy and cost in transient simulations. Parallelization in SEM relies on non-overlapping domain decomposition, partitioning the mesh into subdomains of multiple elements assigned to processors, with inter-processor communication handled via MPI at subdomain interfaces for exchanging GLL point values. This approach achieves excellent weak scalability through localized tensor-product operations within elements, minimizing communication volume to O(P^{d-1}) per interface, and has demonstrated strong scaling efficiency above 90% up to thousands of cores on petascale systems.

Error Analysis

A-priori Error Estimates

In the spectral element method (SEM), a-priori error estimates provide theoretical upper bounds on the approximation error between the exact solution u and its numerical counterpart u_N, where N denotes the polynomial degree within each element. For solutions belonging to Sobolev spaces, the method inherits the standard finite element error bounds, adapted to high-order polynomials. Specifically, for u \in H^{s+1}(\Omega) with $0 \leq s \leq N, the H^1-norm error satisfies

\|u - u_N\|_{H^1(\Omega)} \leq C_s N^{-s} \|u\|_{H^{s+1}(\Omega)},

where C_s is a constant independent of N but dependent on s and the domain \Omega. This algebraic convergence rate highlights the method's ability to achieve high accuracy by increasing N, outperforming low-order methods for sufficiently regular solutions. For analytic solutions, SEM exhibits superior performance due to the spectral nature of the polynomial approximation within elements. The error bound becomes exponential, given by

\|u - u_N\|_{H^1(\Omega)} \leq C \exp(-\gamma N),

where C and \gamma > 0 depend on the analyticity properties of u in a suitable Bernstein ellipse or complex extension of the domain. This near-exponential decay arises from the rapid convergence of spectral projections onto polynomial spaces for smooth, non-periodic functions. These estimates are derived from the Galerkin orthogonality condition inherent in the weak formulation, combined with local approximation theory for polynomial spaces in each element. The orthogonality implies that the error is orthogonal to the finite-dimensional subspace V_N^h, allowing application of the Céa lemma to bound the error by the best approximation error. Local interpolation estimates in Sobolev norms are then extended globally using trace and inverse inequalities to control inter-element continuity, yielding the overall bounds.^[25] The convergence rate depends strongly on the regularity of u: optimal exponential rates are attained for analytic functions, while lower regularity leads to algebraic rates akin to those in finite element methods, though SEM achieves higher orders due to the elevated polynomial degree. For instance, in problems with u \in H^{k+1}(\Omega) for integer k \leq N, the rate improves to O(N^{-k}), surpassing traditional linear elements. Numerical benchmarks confirm these theoretical predictions, particularly the near-exponential convergence for smooth problems. In simulations of 2D elastic wave propagation with analytic initial conditions, SEM solutions demonstrate error reduction by factors of 10 per increment in N by 2-3, aligning closely with the exponential bound and validating the method's efficiency for high-accuracy applications.

Convergence and Stability

The spectral element method (SEM) exhibits high-order convergence properties that depend on the regularity of the solution and the refinement strategy employed. For solutions in the Sobolev space H^{s+1}, the spatial error converges algebraically at a rate of O(N^{-s}), where N is the polynomial degree within each element, assuming sufficient regularity s \geq 1.^[26] For analytic solutions, this improves to exponential convergence, with the error decaying as O(e^{-c N}) under p-refinement (fixed element size, increasing N), enabling rapid accuracy gains even with coarse meshes.^[27] Global convergence follows from aggregating element-wise local estimates, provided the mesh is quasi-uniform to avoid distortion effects that could degrade the rates.^[28] In time-dependent problems, stability for explicit time-stepping schemes is governed by a Courant-Friedrichs-Lewy (CFL) condition, typically requiring \Delta t \leq C \frac{h}{N^2} in one dimension for advection-dominated hyperbolic equations, where h is the element size and C is a constant depending on the wave speed.^[16] This restriction is tighter than in low-order finite element methods due to the high-order nature of SEM, as the minimal inter-node spacing scales as O(h/N^2); however, the use of a diagonal mass matrix from Gauss-Lobatto-Legendre quadrature mitigates this by simplifying inversions and allowing larger effective time steps in practice.^[29] Energy-based stability analyses confirm boundedness in the L^2-norm for hyperbolic systems, with estimates showing \frac{d}{dt} \|u_h\|^2_{L^2} \leq 0 for homogeneous problems on affine meshes, ensuring long-term stability without artificial dissipation.^[28] The h-p refinement strategy in SEM balances element size h and polynomial degree p = N to optimize efficiency, yielding convergence rates such as O(h^{k+1} p^{-m}) for solutions of regularity H^{k+1}, where m relates to the problem's smoothness.^[28] This mixed approach allows exponential overall convergence by increasing p in regions of high regularity while refining h elsewhere, outperforming pure h- or p-refinement for non-uniform solutions. In wave equations, dispersion errors remain minimal due to the orthogonal basis, but high-frequency pollution—accumulative phase errors—can arise; high p (e.g., p \geq 10) significantly reduces this effect, maintaining accuracy with fewer degrees of freedom compared to low-order methods.^[30] Practical enforcement of the CFL condition is essential, often requiring adaptive time-stepping to handle varying N across the domain without compromising stability.^[16]

Historical Development

Origins in the 1980s

The spectral element method (SEM) originated in the mid-1980s as an innovative approach to computational fluid dynamics (CFD), addressing the shortcomings of low-order finite element methods (FEM) in achieving high accuracy for complex flows while overcoming the geometric limitations and ill-conditioning of global spectral methods on non-simple domains. In 1984, Anthony T. Patera introduced the method in his foundational paper, proposing a discretization that partitions the domain into spectral elements—local subdomains where solutions are approximated by high-order polynomials—to blend the versatility of FEM meshing with spectral accuracy. This development was motivated by the need for exponentially convergent solutions to the incompressible Navier-Stokes equations in irregular geometries, such as those encountered in engineering flows, where traditional methods suffered from excessive numerical diffusion or instability.^[31] Patera's early formulation employed a hybrid scheme combining collocation for hyperbolic terms and Galerkin projection for elliptic and pressure terms, using Chebyshev-Gauss-Lobatto points for nodal interpolation within each element to ensure efficient computation and stability. The roots of this approach trace to global spectral methods, which offer superior resolution through high-degree expansions, and local FEM, which handle arbitrary boundaries via element connectivity, thereby resolving spectral methods' sensitivity to domain irregularities. Initial demonstrations focused on two-dimensional incompressible laminar flows, exemplified by simulations of separated flow in a channel expansion, where SEM achieved spectral convergence rates—exponential error reduction with polynomial degree—significantly outperforming standard low-order FEM in accuracy for comparable degrees of freedom.^[31] Key milestones in the late 1980s advanced the method's robustness and applicability. In 1989, Patera and Yvon Maday refined the framework by adopting Legendre polynomials paired with Gauss-Lobatto-Legendre (GLL) quadrature points, which improved mass matrix diagonalization for exact integration and enhanced numerical conditioning over the original Chebyshev basis. Extensions to three-dimensional problems were pursued by Patera and collaborators, enabling tensor-product formulations for 3D Stokes and Navier-Stokes equations while preserving computational efficiency through domain decomposition. Concurrently, integrations with fast iterative solvers, such as preconditioned conjugate gradient methods leveraging sum-factorization techniques, were developed to manage the growing demands of 3D simulations, marking SEM's transition toward practical large-scale CFD implementations.^[32]^[10]

Key Contributors and Advancements

In the late 1980s, Yvon Maday and collaborators refined the Gauss-Lobatto-Legendre (GLL) nodal basis for the spectral element method, demonstrating its advantages in achieving a diagonal mass matrix through Lagrange interpolation and GLL quadrature, which significantly simplifies matrix assembly and inversion in time-dependent simulations.^[8] This work, detailed in their 1989 publication on spectral element methods for the incompressible Navier-Stokes equations, established key theoretical properties that enhanced the method's efficiency for high-order approximations.^[33] During the 1990s, Anthony T. Patera and his group at MIT advanced the spectral element method through developments in hp-refinement strategies, which combine h-refinement (mesh size reduction) with p-enrichment (polynomial degree increase) to optimize accuracy and computational cost, and robust preconditioners to accelerate iterative solvers for large-scale problems. In the late 1990s, Dimitri Komatitsch and Jeroen Tromp extended the method to three-dimensional seismic wave propagation, facilitating simulations in realistic Earth models.^[34] These contributions culminated in the creation of the Nekton and subsequent Nek5000 codes, open-source frameworks that implemented parallel spectral element solvers for fluid dynamics and beyond, enabling efficient simulations on distributed-memory architectures.^[35] A key milestone of this era was the focus on parallel implementations, such as tuned spectral element methodologies that facilitated load balancing and domain decomposition on early parallel machines, achieving scalable performance for elliptic and time-dependent problems.^[36] In the 2000s, Jan S. Hesthaven and colleagues extended the spectral element framework to discontinuous Galerkin spectral element methods (DGSEM), particularly for hyperbolic problems, by incorporating nodal collocation on GLL points within discontinuous subspaces to handle discontinuities and improve stability without global continuity enforcement. This advancement, as outlined in their comprehensive 2008 monograph on nodal DG methods, provided high-order accuracy and energy stability for applications involving wave propagation and advection-dominated flows.^[37] Post-2010 developments have emphasized computational efficiency and adaptability, with Dimitri Komatitsch leading efforts in GPU acceleration for seismic wave propagation using spectral elements, achieving speedups of 20 to 60 times over CPU implementations through hybrid CPU-GPU algorithms that preserve the method's high accuracy.^[38] His work, including fluid-solid coupling simulations on multi-GPU clusters starting around 2010, has enabled large-scale modeling of complex geological structures.^[38] Concurrently, adaptive techniques have evolved, building on earlier hp-refinements to dynamically adjust polynomial orders and mesh resolutions based on error estimators, further integrated in codes like Nek5000 for real-time optimization in multiphysics simulations.^[39]

Applications

Computational Fluid Dynamics

The spectral element method (SEM) has been extensively applied to the incompressible Navier-Stokes equations, leveraging its high-order accuracy to simulate transitional and turbulent flows with superior resolution compared to low-order methods. This approach excels in capturing complex vortical structures and boundary layers in viscous flows, as demonstrated in simulations of pipe flows using the open-source solver Nek5000, where SEM enables direct numerical simulations (DNS) at Reynolds numbers up to 10^4 with polynomial orders of 7-15, achieving exponential convergence within each element.^[40]^[41] To enforce the incompressibility constraint in SEM discretizations, adaptations such as pressure projection methods are commonly employed, where an intermediate velocity field is projected onto a divergence-free space via a Poisson equation for pressure correction, ensuring mass conservation without introducing spurious modes.^[42] Alternatively, divergence-free bases, constructed using curl-conforming elements or specific polynomial spaces, can directly parameterize solenoidal velocity fields, reducing the need for iterative pressure solves in certain formulations.^[33] For compressible flows with variable density, SEM extends to the full Navier-Stokes equations by incorporating entropy-stable flux formulations and handling density variations through equation-of-state models, as seen in high-fidelity simulations of reacting flows.^[43] Benchmark problems like the Taylor-Green vortex decay illustrate SEM's spectral-like convergence, where energy dissipation rates match reference DNS solutions with errors below 0.1% using modest grids of 8^3 elements at order 15, highlighting its efficiency for transitional regimes. In large-eddy simulations (LES) of turbulent flows, SEM requires fewer elements than traditional DNS—often 10-100 times less—due to its high-order resolution of large scales, enabling under-resolved modeling of subgrid stresses via explicit filtering.^[44]^[45] The method's efficiency stems from reduced degrees of freedom (DOFs), allowing high-Reynolds-number simulations on parallel clusters; for instance, Nek5000 has enabled petascale simulations of nuclear fuel assemblies and coolant flows at Re up to ~80,000, achieving high performance with matrix-free operators. Its GPU-accelerated successor, NekRS, further enables exascale simulations of complex flows, including nuclear reactor components on systems like Frontier.^[46]^[47]^[48] However, challenges arise at inflow and outflow boundaries, where non-reflective conditions are essential to prevent spurious reflections; special treatments like characteristic projections decompose waves into acoustic, vorticity, and entropy modes, projecting only outgoing components to maintain stability in open-domain flows.^[49]

Seismic and Wave Propagation

The spectral element method (SEM) is extensively applied in seismic modeling to solve the 3D elastic wave equation in complex, heterogeneous global Earth models, enabling accurate simulations of earthquake wave propagation. The open-source SPECFEM3D software package implements SEM for such computations, facilitating forward modeling of seismic events from source mechanisms to surface ground motions. This approach is particularly valuable for global-scale simulations, where it handles irregular geometries and material discontinuities inherent in Earth's crust, mantle, and core.^[50] A key advantage of SEM in seismic applications is its low numerical dispersion, which preserves the accuracy of high-frequency waves over long propagation distances, outperforming lower-order methods in resolving fine-scale features. Additionally, SEM provides precise implementation of free-surface boundary conditions and efficient absorbing boundaries through perfectly matched layers (PML), minimizing spurious reflections in unbounded domains. These properties make SEM ideal for hyperbolic wave problems in viscoelastic media, with high-order polynomial bases (degree p = 8-15) achieving resolutions up to 5 Hz for crustal-scale studies. In practice, SEM has been used to model tsunami propagation by coupling seismic sources with ocean acoustics, as seen in simulations of tsunamigenic earthquakes that track wave energy transfer from solid Earth to water columns. For crustal imaging, SEM-based full-waveform inversion (FWI) reconstructs subsurface velocity structures from observed seismograms, enhancing resolution in fault zones and sedimentary basins.^[51] Adaptations include the use of curvilinear hexahedral elements to conform to spherical domains, ensuring geometric fidelity in global models, and integration with frequency-dependent attenuation models via generalized standard linear solids to account for anelastic losses. Recent advancements in the 2020s leverage SEM for real-time earthquake forecasting, such as the Real-time Online earthquake Simulation (ROS) system in Taiwan, which uses SPECFEM3D to generate ground-motion predictions within minutes of event detection, supporting rapid hazard assessment.^[52] Exascale computing has enabled large-scale FWI applications, with simulations on supercomputers like ORNL's Frontier processing thousands of events to achieve mantle-scale imaging at periods down to 17 seconds, paving the way for higher-frequency global inversions.^[53]^[54]

Spectral Methods

Global spectral methods approximate solutions to partial differential equations using basis functions, such as Fourier series for periodic domains or Chebyshev polynomials for non-periodic but simple geometries, that extend across the entire computational domain. These methods excel in achieving high accuracy for smooth solutions on regular domains, leveraging the global nature of the basis to attain exponential convergence rates. However, they are susceptible to aliasing artifacts, particularly in nonlinear problems, which necessitate techniques like dealiasing filters or over-integration to maintain stability. The spectral element method (SEM) differs by decomposing the domain into a mesh of non-overlapping elements, where high-order polynomial expansions—typically tensor-product bases—are applied locally within each element via mappings to reference domains. This local approximation enables SEM to accommodate complex geometries and irregular boundaries that challenge global spectral methods, though it introduces inter-element continuity enforcement through coupling at interfaces. Despite this added complexity, SEM preserves the spectral accuracy of exponential convergence within individual elements for smooth functions.^[31] A key trade-off lies in computational efficiency: global spectral methods are more cost-effective for smooth, periodic problems on simple domains, scaling as O(N^d) operations per time step in d dimensions, where N is the polynomial degree, due to fast transform algorithms like the FFT. In SEM, the cost rises to approximately O(E N^{d+1}), with E the number of elements, owing to the overhead of local operations and global assembly, making it less efficient for straightforward cases but more versatile overall.^[21] SEM serves as a bridge between global spectral methods and traditional finite elements; when the domain is covered by a single element, it reduces exactly to a global spectral approximation. In practice, SEM is often preferred over pseudospectral (global) methods for non-periodic computational fluid dynamics simulations involving irregular boundaries, such as flows in channels with expansions or around obstacles, where global bases fail to adapt without excessive distortion.^[31]

Finite Element Methods

The standard finite element method (FEM) employs low-order polynomial basis functions, typically linear (p=1) or quadratic (p=2) elements, to approximate solutions over a domain divided into a mesh of simple geometric shapes such as triangles or tetrahedra in two or three dimensions. Accuracy in FEM is primarily achieved through h-refinement, which involves reducing the size of mesh elements (decreasing h) to increase resolution, leading to algebraic convergence rates of order O(h^{p+1}) for smooth solutions. This approach offers versatility for arbitrary, unstructured meshes, making it suitable for complex geometries where mesh conformity to boundaries is essential.^[55] In contrast, the spectral element method (SEM) emphasizes high-order p-refinement, using polynomials of degree p=4 to 15 or higher within larger, fixed-size elements, which enables exponential convergence for smooth problems while keeping the coarse mesh resolution (h) constant. This high p-order in SEM yields superior efficiency for problems requiring high accuracy, as fewer degrees of freedom are needed compared to the algebraic scaling in low-order FEM.^[55] However, SEM relies on tensor-product elements, typically quadrilaterals or hexahedra, which restricts its application to more structured or semi-structured domains, unlike the simplex-based flexibility of standard FEM. Both SEM and FEM are formulated within the Galerkin framework, projecting the governing equations onto finite-dimensional subspaces, but they differ in basis representation and matrix properties.^[55] SEM uses nodal Gauss-Lobatto-Legendre (GLL) basis functions for efficient collocation, whereas standard FEM employs Lagrange or hierarchical bases on simplex elements. A key distinction lies in the mass matrix: SEM achieves a diagonal mass matrix through exact GLL quadrature, facilitating explicit time-stepping, while FEM often uses lumped (diagonalized approximation) mass matrices for similar purposes but with less exactness.^[55] FEM is generally preferred for problems involving irregular geometries or solutions with low regularity, where h-refinement allows adaptive meshing around singularities without high polynomial costs. Conversely, SEM excels in scenarios demanding high accuracy for smooth fields, such as in computational fluid dynamics with periodic or structured domains, due to its rapid p-convergence.^[55] SEM can be viewed as a specialized case of the more general hp-FEM, where uniform high-p refinement is applied across tensor-product elements, combining the exponential accuracy of spectral methods with the domain decomposition of finite elements. This positions SEM within the broader hp framework, which allows variable p and h adaptations but simplifies to fixed h and high uniform p in the spectral variant.^[55]

References

[1]
[PDF] Introduction to Finite and Spectral Element Methods using Matlab
What is the spectral element method? The spectral element method is an advanced implementation of the finite ele- ment method in which the solution over ...
[2]
[PDF] The Spectral-Element Method in Seismology
The spectral-element method is for 3D Earth models, using hexahedral elements, combining finite-element and pseudospectral methods, and using a weak ...
[3]
https://arxiv.org/pdf/1802.06743.pdf
[4]
[PDF] A Review: Applications of the Spectral Finite Element Method
Jan 11, 2023 · Moreover, Fig. 5 shows the applica- tions of the spectral element method within specific research areas. In Sect. 3, work on various strategies ...
[5]
A spectral element method for fluid dynamics: Laminar flow in a ...
A spectral element method that combines the generality of the finite element method with the accuracy of spectral techniques is proposed
[6]
Spectral Element Method - an overview | ScienceDirect Topics
The spectral element method is defined as a Galerkin method that utilizes a finite dimensional space of polynomials, whose degree does not exceed m in each ...
[7]
[PDF] The Spectral-Element Method
Originated in fluid dynamics (Patera, 1984, Maday and Patera, 1989). • First applications to seismic wave propagation by Priolo, Carcione and Seriani, 1994).
[8]
[PDF] THE SPECTRAL ELEMENT METHOD (SEM): FORMULATION AND ...
Patera, A.T. (1984). A Spectral element method for fluid dynamics: laminar flow in a channel expansion, Journal of. Computational Physics, 54, 468-488.
[9]
A Review: Applications of the Spectral Finite Element Method
Mar 28, 2023 · The benefit of using spectral elements is that they may be used in various scenarios to produce steady computational solutions and high accuracy ...
[10]
[PDF] Spectral Element Methods: Algorithms and Architectures
In the spectral element discretization (Patera, 1984; Maday and Patera, 19871, the computational domain i s broken up into macro-spectral elements, and the ...
[11]
[PDF] Efficiency of the spectral element method with very high ... - HAL
The Spectral Element Method (SEM) has gained tremendous popularity within the seismological community to solve the wave equation at all scales.
[12]
Spectral Element Methods for Elliptic Problems in Nonsmooth ...
The accuracy of the spectral element methods is then degraded, and they offer no apparent advantage over low-order finite element methods. In many cases ...Missing: limitations | Show results with:limitations
[13]
https://arxiv.org/pdf/1706.00855.pdf
[14]
MASSIVELY-PARALLEL SPECTRAL ELEMENT ... - OAKTrust
tensor product-based quadrilateral and hexahedral spectral elements to more general shapes, such as triangle and tetrahedrals, and give examples for Euler ...
[15]
An efficient implicit spectral element method for time-dependent ...
We present an implicit spectral element method to approximate the solution of time-dependent nonlinear diffusion equations in complex geometries.
[16]
[PDF] Stability of an explicit high-order spectral element method for ... - HAL
Oct 18, 2018 · spectral element method ... matrices are computed in the reference element after introducing the isoparametric mapping and a numerical quadrature.
[17]
[PDF] Contents
In the spectral element method (SEM), data is represented on sets of non ... are the coordinates in the canonical reference element, Ω := [−1,1]3.
[18]
[PDF] Spectral Equivalence Properties of Higher-Order Tensor ... - OSTI
• Tensor product elements (quadrilateral & hexahedral). • Span same polynomial spaces as other “standard” shape functions. • Enable correspondence between ...
[19]
[PDF] adaptive mesh strategies for the spectral element method
In this paper, we present an adaptive spectral element method which automatically allocates resolution where it is most needed in an optimal fashion.
[20]
Adaptive mesh strategies for the spectral element method
Adaptive mesh strategies that include resolution refinement and coarsening by three different methods are illustrated on solutions to the one-dimensional ...
[21]
[PDF] From h to p Efficiently: Implementing finite and spectral/hp element ...
The spectral/hp element method can be considered as bridging the gap between the. – traditionally low order – finite element method on one side and spectral ...
[22]
[PDF] Legendre Spectral Finite Elements for Reissner-Mindlin Plates
This choice of element nodes has two important advantages. First, the uneven spacing of the GLL points causes a “clustering” of element nodes at element ...
[23]
A two-dimensional spectral-element method for computing spherical ...
For clarity, we shall denote GLL nodes in the ~ direction as ~ and in the i~ direction as i~ . The choice of Lagrangian interpolation upon GLL nodes yields ...
[24]
Stability of an explicit high‐order spectral element method for ...
Jul 4, 2018 · ... spectral element method ... The elemental matrices are computed in the reference element after introducing the isoparametric mapping and a ...
[25]
A Class of Spectral Element Methods and Its A Priori/A Posteriori ...
Oct 27, 2013 · The special work of this paper is as follows. (1) We prove a priori and a posteriori error estimates for spectral and spectral element methods.Missing: seminal | Show results with:seminal
[26]
[PDF] Stability estimates for hp spectral element methods
Abstract. In this paper we show that the h-p spectral element method developed in [3,8,9] applies to elliptic problems in curvilinear polygons with mixed ...
[27]
[PDF] 2D Spectral Element Scheme for Viscous Burgers' Equation
May 18, 2004 · To solve these equations efficiently while maintaining a high working accuracy over long time periods, we choose the Spectral Element Method ( ...<|control11|><|separator|>
[28]
[PDF] an analysis of a non-conforming hp-discontinuous galerkin spectral ...
We analyze the consistency, stability, and convergence of a hp-discontinuous Galerkin spectral element method. The analysis is simultaneously done for ...
[29]
[PDF] Notes on the implementation of spectral element method (SEM) for ...
Spectral element method (SEM) is a specific form of finite element method (FEM), where the mass matrix can be trivally inverted thanks to its diagonal ...
[30]
[PDF] Efficient implementation of high-order finite elements for Helmholtz ...
... pollution effect and as a result their accuracy deteriorates as the frequency ... Finite element solution to the Helmholtz equation with high wave number – Part ...
[31]
https://doi.org/10.1016/0021-9991(84)90128-1
[32]
Spectral element methods for the incompressible Navier-Stokes ...
The theoretical foundations and numerical implementation of spectral element methods for the incompressible Navier-Stokes equations are presented, considering ...
[33]
Spectral element methods for the incompressible Navier-Stokes ...
Spectral element methods for the incompressible Navier-Stokes equations · Y. Maday, A. Patera · Published 1989 · Engineering, Physics.
[34]
[PDF] Spectral/hp Element Framework - Nektar++
Jan 20, 2021 · 2.2.2.5 The spectral element method. Patera [34] combined the high accuracy of the spectral methods with the geometric flexibility of the ...
[35]
[PDF] Nek5000 Tutorial - Mathematics and Computer Science
▫Parallelism and geometric flexibility derive from the unstructured, element-by-element, operator evaluation. ▫Elements, or groups of elements are distributed ...Missing: parallelization | Show results with:parallelization
[36]
A spectral element methodology tuned to parallel implementations
A spectral element method for fluid dynamics: laminar flow in a channel expansion. J. Comput. Phys. (1984) ; Parallel spectral element solution of the Stokes ...
[37]
Fluid–solid coupling on a cluster of GPU graphics cards for seismic ...
We develop a hybrid multiGPUs and CPUs version of an algorithm to model seismic wave propagation based on the spectral-element method in the case of models ...
[38]
[PDF] Aspects of adaptive mesh refinement in the spectral element method
Jun 11, 2019 · It is also shown that, on a single core, Nek5000 is memory limited rather than compute limited. Then, a new method for the coarse grid part of ...
[39]
A Review on Spectral Element Solver Nek5000 - AIP Publishing
Basically Nek5000 based on spectral elements code, Nekton 2.0. For distributed-memory parallel processors Nekton 2.0 was the first available three dimensional ...Missing: parallelization | Show results with:parallelization
[40]
[PDF] Turbulence Modeling with Nek5000/RS, SOD2D and Alya
Sep 25, 2023 · We demonstrate the accuracy of these codes for two benchmark sub-sonic turbulent flows: periodic channel and pipe flow, by carrying out first ...
[41]
[PDF] A spectral element projection scheme for incompressible flow with ...
A projection scheme based on the pressure correction method is discussed to solve the Navier-Stokes equations for incompressible flow.
[42]
A discontinuous Galerkin spectral element method for compressible ...
In this work, a high-order discontinuous Galerkin spectral element method (DGSEM) is developed to solve the chemically reacting Navier-Stokes equations.
[43]
[PDF] Simulation of the Taylor–Green Vortex Using High-Order Flux ...
The reference data (solid black contours) were computed using a spectral element method with 5123 DOFs by Carton de Wiart et al. [37] and are available from the ...
[44]
[PDF] Spectral element filtering techniques for large eddy simulation with ...
Spectral element methods have previously been successfully applied to direct numerical simulation of turbulent flows with moderate geometrical complexity ...
[45]
Energy Exascale Computational Fluid Dynamics Simulations With ...
Spatial discretization is based on the high-order spectral element method (SEM), which affords the use of fast, low-memory, matrix-free operator evaluation.
[46]
Spectral element applications in complex nuclear reactor geometries
The spectral element code Nek5000 is an open-source, higher-order computational fluid dynamics code developed at Argonne National Laboratory.
[47]
[PDF] Strong compact formalism for characteristic boundary conditions ...
Feb 24, 2020 · Overall, this work enables the use of strong discontinuous spectral methods to study unsteady problems on complex geometries. Keywords: Spectral ...<|control11|><|separator|>
[48]
SPECFEM | A family of open-source spectral-element method ...
The SPECFEM codes use the spectral-element method to simulate seismic wave propagation on different scales. Together with dedicated inversion tools, they ...Missing: real- time forecasting 2020s exascale full-
[49]
Spectral Element Modeling of Acoustic to Seismic Coupling Over ...
Dec 22, 2021 · We use the SPECFEM3D software to model the propagation and coupling of acoustic and seismic waves with the spectral element method (Komatitsch ...
[50]
Toward real-time regional earthquake simulation II
Jun 15, 2014 · ... spectral-element method (SEM) to simulate seismic wave propagation. The SEM is a numerical technique developed 30 years ago to address ...
[51]
[PDF] Global Full-Waveform Inversion: Towards Exascale Imaging of ...
hexahedral elements". Gauss-Lobatto-Legendre quadrature" finite elements with high-degree polynomial interpolation: accuracy of pseudo-spectral.
[52]
Spectral/hp Element Methods for Computational Fluid Dynamics
Free delivery 25-day returnsThis text aims to introduce a wider audience to the use of spectral/hp element methods with particular emphasis on their application to unstructured meshes.