PaperlayerView exact revision on arXiv

Source-linked AI summary

$p$ and $rp$-Spectral Element Methods for Two-dimensional Elliptic Boundary Layer Problems

Sonia, Akhlaq Husain, Arshad Khan

arXiv:2609.29185v1math.NAcs.NA

TL;DR

Elliptic boundary-layer problems develop sharp gradients near boundaries, challenging standard numerical methods. The paper proposes p- and rp-least-squares spectral element methods for rectangular domains, using fixed or layer-adapted meshes. The rp-version has a convergence rate uniform in ε, and numerical results support the methods’ ability to resolve layers.

  • Problem

    Elliptic boundary-layer problems involve sharp gradients in narrow regions near the boundary, which standard methods can struggle to resolve accurately.

    Source: §1 · pp. 1–2

  • Method

    The paper develops non-conforming p- and rp-least-squares spectral element methods, with stability estimates and preconditioning for the resulting normal equations.

    Source: §1 · p. 2, §5 · p. 11, §8 · p. 28

  • Results

    The rp-version achieves a rate O(εα^(2W)) uniformly in ε, and numerical results confirm the method’s efficiency and ability to resolve boundary layers on rectangular domains.

    Source: §8 · p. 29

  • Takeaways & Limitations

    The numerical results validate the theoretical error estimates and predicted computational complexity, and the method is practically implementable on parallel computers.

    Source: §8 · p. 29

    The stated exponential convergence for the layer-focused mesh requires smooth solution components to be approximated at an exponential rate as well.

    Source: §6 · p. 16

Abstract

from arXiv · show

We propose $p$ and $rp$-spectral element methods for elliptic boundary layer problems on two dimensional rectangular domains. To resolve boundary layers, we use $p$-version with a fixed mesh and an $rp$-version with a boundary layer mesh consisting of thin needle-like elements near the boundary layer and coarse elements away from the layer. Stability estimates are derived using non-conforming spectral element functions. The numerical scheme for both the methods is based on minimising a least-squares functional in appropriate Sobolev norms. We construct robust preconditioners to manage the condition number of the normal equations derived from the least-squares formulation. The method is able to approximate boundary layers at a rate $O\left(\frac{\log W}{W^2}\right)$ for the $p$-version and at the rate $O({\epsilon}\alpha^{2W})$, uniformly in $\epsilon$, for the $rp$-version, where $0<\epsilon\leq1$ is the boundary layer parameter, $\alpha<1$ is a constant and $W$ denotes the degree of the approximating polynomial. Numerical results are provided for model elliptic boundary layer problems for a range of boundary layer parameters. Simulation results demonstrate the efficiency and robustness of the method in capturing boundary layers on rectangular domains.

1. Introduction

Elliptic boundary layers create sharp gradients near domain boundaries, which can challenge standard numerical methods. The paper proposes p- and rp-least-squares spectral element methods for rectangular domains, with the rp-version using a boundary-layer mesh.

1. Introduction

  • Boundary layers are narrow regions of sharp gradients, and standard numerical methods may need boundary refinement to resolve them accurately.The introduction notes that some methods may not converge in rare scenarios.Source: §1 · pp. 1–2
  • Prior approaches include h-, p-, and hp-refinement, layer-adapted meshes, and high-order spectral element methods.The paper motivates least-squares spectral elements by combining the accuracy of spectral methods with the flexibility of finite elements.Source: §1 · p. 2
  • The proposed methods minimize least-squares residuals in fractional Sobolev norms and penalize jumps across element interfaces to enforce continuity.They target elliptic boundary-layer problems on two-dimensional rectangular domains with Dirichlet boundary conditions.Source: §1 · p. 2
  • The rp-version achieves a convergence rate O(εα^W), uniformly in ε, on a boundary-layer mesh.Here 0 < ε ≤ 1 is the boundary-layer parameter and α < 1 is a constant.Source: §1 · p. 2

2. Elliptic boundary layer problem

The section introduces Sobolev-space notation and a rectangular elliptic boundary-layer model. Its analysis focuses on the dominant exponential decay normal to the boundary, while the solution may also contain smooth and corner-layer components.

2.1. Notations and preliminaries

  • The notation section introduces integer-order Sobolev spaces, zero-trace functions, fractional Sobolev norms, and a parameter-weighted norm.It also establishes notation for points and derivatives in two dimensions.Source: §2.1 · p. 3

2.2. The model problem

  • The model is posed on the square Ω = (−1, 1)² with Dirichlet boundary conditions.The square has a Lipschitz boundary and corners.Source: §2.2 · p. 4
  • The solution can contain regular layers along the four smooth edges and corner layers at the vertices, alongside smooth components.The regularity estimate is non-uniform in ε because its constant depends on ε.Source: §2.2 · pp. 4–5
  • The theoretical analysis targets the dominant one-dimensional exponential decay in the direction normal to the boundary.The layer profile is considered away from corners, with ρ denoting normal distance and θ the coordinate along an edge.Source: §2.2 · p. 4

3. Domain Discretization and stability estimate

The domain is partitioned into quadrilateral elements represented by tensor-product polynomials mapped from a reference square. The stability estimate has an ε-dependent constant for the fixed-mesh p-version, motivating the rp-version’s boundary-layer mesh.

3.1. Discretization and spectral element representation

  • Each quadrilateral element uses tensor-product polynomials mapped from the bi-unit square to the physical domain.The discretization allows uniform or non-uniform refinement.Source: §3.1 · pp. 5–6

3.2. Stability estimate

  • Jumps in the function and its derivatives are defined across shared element edges as part of the stability framework.The passages describe the shared-edge mappings and derivative jumps.Source: §3.2 · p. 6
  • The stability estimate’s constant depends on ε but not W, making the fixed-mesh p-version estimate non-uniform in ε.The text identifies this limitation as motivation for the rp-version.Source: §3.2 · pp. 7–8

4. Numerical scheme

The numerical scheme minimizes a least-squares functional that combines equation and boundary residuals with inter-element jumps. The resulting normal equations are solved with PCGM.

4. Numerical scheme

  • The discrete solution minimizes the functional over the spectral element function space, and PCGM solves the resulting normal equations.The normal equations are written as AT AU = AT b.Source: §4 · p. 8
  • The least-squares functional combines PDE residuals, boundary-condition residuals, and jumps in the function and its derivatives across element interfaces.Boundary residuals and jumps are measured in appropriate fractional Sobolev norms; jump terms enforce inter-element continuity.Source: §4 · p. 8
  • PCGM iterations require matrix-vector products, which can be computed without explicitly storing mass and stiffness matrices.Neighboring elements exchange function values and derivatives at interfaces during iterations.Source: §4 · p. 9

5. Preconditioning techniques and parallelization

The method uses a scaled H2 quadratic form as a preconditioner for the least-squares system. Its elementwise structure supports efficient inversion and yields a condition number of O((log W)^2).

5. Preconditioning techniques and parallelization

  • The preconditioner is spectrally equivalent to the least-squares quadratic form, with equivalence constants independent of N and W but dependent on ε.The text describes the preconditioner as computationally easy to invert and designed for robustness across parameters.Source: §5 · p. 9
  • The elementwise scaled H2 quadratic form preconditions the least-squares form, with a resulting condition number of O((log W)^2).The preconditioner is composed of decoupled quadratic forms on each element.Source: §5 · pp. 9, 11
  • A polynomial basis diagonalizes the elementwise preconditioner, and PCGM communication is limited to interface values and derivatives plus two global scalars.The text reports approximately O(W ln W) PCGM iterations for exponential accuracy.Source: §5 · p. 11

6. Error estimates

The error estimates distinguish a fixed-mesh p-version from an rp-version using a boundary-layer mesh. The rp-version gives uniform exponential convergence for the boundary-layer component, subject to resolution of any smooth components.

6. Error estimates

  • The fixed-mesh p-version requires W = O(1/ε) to resolve the layer, so fixed W leads to degraded accuracy as ε decreases.The text describes the strict asymptotic requirement as inherently non-uniform for the p-version.Source: §6 · p. 15
  • The two-element mesh with O(1) polynomial degree away from the layer may not approximate smooth solution components exponentially; those components need higher degree or further subdivision.Overall exponential convergence requires the smooth components to be approximated at an exponential rate.Source: §6 · pp. 16–17
  • A two-element boundary-layer mesh gives exponential convergence uniformly in ε for the boundary-layer problem.The mesh uses a thin strip near the layer and a coarse element away from the boundary.Source: §6 · pp. 15–16

7. Numerical results

Numerical experiments compare p-refinement on uniform meshes with rp-refinement on meshes aligned to boundary layers, across several layer locations and values of ε. Uniform p-refinement improves with mesh refinement but remains limited for the sharpest layers, while the layer-aligned rp examples show rapid error reduction.

Numerical setup

  • The Example 7.1 p-version tests used uniform meshes with 1, 2, 4, and 9 elements, degrees W=2, 4, …, 32, and ε values from 1 to 0.01.The section describes simulations using FORTRAN-90 and a 120-core system.Source: §7 · pp. 17–19

Example 7.1

  • For Example 7.1, uniform p-refinement improved with finer meshes, but convergence remained suboptimal for ε≤0.025, motivating mesh adaptation or rp-refinement.On a single element, errors saturated or increased for small ε; a nine-element mesh improved results, especially for moderate ε.Source: §7 · p. 19

Example 7.2

  • In Example 7.2, four elements yielded exponential convergence for moderate and small ε, while nine elements achieved robust exponential convergence even at ε=0.01.The single-element mesh stagnated for smaller ε, and two elements remained inadequate for strong layers.Source: §7 · pp. 21–22

Example 7.3

  • For Example 7.3, nine elements gave robust exponential convergence for ε≥0.075 and reduced error even at ε=0.01, although sharper layers still call for layer-adapted methods.A single-element mesh failed to resolve small-ε layers, while two- and four-element meshes captured moderate layers more effectively.Source: §7 · pp. 23–25

Examples 7.4–7.5

  • On layer-aligned rp meshes, for ε = 0.01, error decreases rapidly once W ≥ 20; a reported error value is 5.56 × 10^-2 at W = 32.The authors say these results validate the theoretical estimate in Theorem 6.2.Source: §7 · pp. 26, 28

8. Conclusion and future prospects

The paper presents a non-conforming least-squares spectral element method for two-dimensional elliptic boundary layer problems on rectangular domains, using fixed-mesh p-refinement and layer-adapted rp-refinement. It reports approximation rates for both versions and says numerical results validate the theoretical estimates and predicted computational complexity.

Conclusion and future prospects

  • Both versions use a least-squares functional in appropriate Sobolev norms, with stability estimates based on non-conforming spectral element functions.The method also uses robust preconditioners to control the condition number of the normal equations associated with the least-squares formulation.Source: §8 · pp. 28–29
  • The rp-version achieves an approximation rate of O(epsilon alpha^2W), uniformly in epsilon, where 0 < epsilon <= 1 is the boundary-layer parameter, alpha < 1 is a constant, and W is the approximating-polynomial degree.Source: §8 · p. 29
  • Numerical results validate the theoretical error estimates and predicted computational complexity; the method is also described as practically implementable on parallel computers.The authors identify smooth domains with analytic boundaries as a topic for future work.Source: §8 · p. 29

Funding

The authors report receiving no financial support for the submitted work.

Funding

  • The authors received no financial support from any organization for the submitted work.Source: Funding · p. 29

Data Availability Statement

The article states that it contains all data generated or analyzed in the study.

Data Availability Statement

Declaration of AI use

The authors state that no AI-assisted technologies were used in preparing the article.

Declaration of AI use

Appendix A. Technical Results

Appendix A constructs vertex, side, and edge corrections for piecewise spectral element functions, yielding a globally H2 function while preserving the boundary data. It then uses that construction to establish a further estimate for the corrected function.

Appendix A. Technical Results

  • The correction in Lemma A.1 produces a globally H2 function u + v, with the elementwise corrections zero on Γ and a constant Cǫ independent of W.The excerpt introduces Cǫ as depending on ǫ; the estimate itself is not visible here.Source: Appendix A · pp. 32, 34
  • Vertex corrections use averages over elements sharing each vertex to set polynomial values and first derivatives at the four square vertices.The polynomial is constructed to match the prescribed values and derivatives at each vertex.Source: Appendix A · p. 32
  • Side corrections match averaged function values and first derivatives across elements sharing a side, while vanishing at the square’s vertices.For edges contained in Γ ∪ Γ0, the subsequent edge-correction terms are set identically to zero.Source: Appendix A · p. 33
  • Lemma A.2 applies Lemma A.1 to w = u + v, which belongs to H2(Ω) and agrees with u on Γ, to obtain its stated estimate.The estimate’s expression is not present in the supplied passage.Source: Appendix A · p. 34
Loading 2609.29185v1…