Optimal error estimates for a discontinuous Galerkin method on curved boundaries with polygonal meshes

Technical Summary: Optimal Error Estimates for a Discontinuous Galerkin Method on Curved Boundaries with Polygonal Meshes

Problem Statement
The paper addresses the numerical solution of boundary value problems, specifically two-dimensional linear advection-diffusion-reaction equations, on domains with arbitrary smooth curved boundaries ($\Omega$). A primary challenge in this context is the loss of convergence order that occurs when the physical domain is approximated by a polygonal mesh ($\Omega_h$). Standard Discontinuous Galerkin (DG) methods, when applied to such polygonal approximations with boundary conditions imposed directly on the computational boundary ($\partial\Omega_h$), generally suffer from suboptimal convergence. The geometric error introduced by the mismatch between the curved physical boundary and the piecewise linear computational boundary can pollute the solution, limiting accuracy to at most second-order regardless of the polynomial degree used.

Methodology
The authors propose and analyze the DG–ROD (Discontinuous Galerkin – Reconstruction for Off-site Data) method. This approach extends a strategy previously developed for Finite Element (FE) and Finite Volume (FV) methods to the DG framework.

• Mesh and Space: The physical domain is discretized using standard polygonal (triangular) meshes. The method utilizes broken Sobolev spaces and discontinuous piecewise polynomial spaces ($S_{h,N}$).
• Boundary Condition Treatment: Instead of imposing the prescribed boundary conditions directly on the polygonal boundary $\partial\Omega_h$, the DG–ROD method employs a polynomial reconstruction. For each boundary element, a constrained minimization problem is solved to determine a local polynomial that minimizes the distance to the numerical solution while satisfying the exact boundary conditions at specific points on the true physical boundary $\partial\Omega$.
• Variational Formulation: The method combines the Symmetric Interior Penalty (SIP) scheme for diffusion and an upwind scheme for advection-reaction terms. The reconstructed boundary data ($g_{ROD}$) is incorporated directly into the primal variational formulation, effectively correcting the geometric mismatch within the weak form.
• Implementation: The authors adopt a global formulation where the reconstruction and the DG solver are solved in a single step, rather than an iterative alternating approach.

Key Contributions and Theoretical Analysis
The paper provides a rigorous theoretical analysis of the DG–ROD method, establishing existence, uniqueness, and optimal error estimates for both convex and non-convex domains.

1. Well-Posedness: The authors prove the boundedness of the bilinear form and establish an inf-sup condition (via the Banach–Nečas–Babuška theorem), ensuring the existence and uniqueness of the discrete solution for sufficiently small mesh size $h$ and sufficiently large penalty parameter $\eta$.
2. Convex Domains: For convex domains, the authors prove that the Galerkin orthogonality holds. Under suitable regularity assumptions on the exact solution, they derive:
   • Optimal convergence of order $O(h^N)$ in the DG norm (energy norm).
   • Optimal convergence of order $O(h^{N+1})$ in the $L^2$ norm using duality arguments.
3. Non-Convex Domains: In non-convex settings, the geometric approximation breaks the standard Galerkin orthogonality. The authors introduce a residual term arising from the extension of the solution to a larger smooth domain $\tilde{\Omega}$. They derive error estimates that include a residual term dependent on the geometric mismatch.
   • Theoretical analysis suggests that for non-convex domains with $N > 2$, optimal convergence in the $L^2$ norm may not be achieved unless the computational domain approximates the physical domain with higher accuracy (order $h^q, q>2$) or the solution admits a regular extension satisfying specific conditions.
   • The analysis highlights that the loss of orthogonality constitutes the main technical challenge in the non-convex case.

Numerical Results
The theoretical findings are validated through numerical experiments on three benchmark geometries: a disk (convex), an annulus (non-convex), and a rose-shaped domain (non-convex).

• Comparison: Results compare the classical DG method (imposing conditions on $\partial\Omega_h$) against the DG–ROD method.
• Performance: The classical DG method consistently demonstrates accuracy deterioration, with convergence rates limited to approximately second order, regardless of the polynomial degree $N$.
• DG–ROD Efficacy: The DG–ROD method successfully recovers optimal convergence rates ($O(h^{N+1})$ in $L^2$) for $N=2, 3, 4$ across all tested geometries and various advection-reaction coefficients. Notably, even in the non-convex rose-shaped domain, the method achieves optimal rates, which the authors attribute to the use of smooth manufactured solutions that allow for regular extensions, effectively bypassing the theoretical worst-case limitations derived for general non-convex cases.

Significance and Claims
The paper claims to provide the first rigorous theoretical analysis of the DG–ROD method, extending the theoretical framework from Finite Elements to Discontinuous Galerkin methods.

• Restoration of Accuracy: The primary significance is the demonstration that optimal convergence rates can be restored on simple polygonal meshes without the need for complex curved elements or isoparametric mappings, provided the boundary data is reconstructed via the ROD technique.
• Robustness: The numerical results confirm the method's robustness in handling both convex and non-convex domains with curved boundaries.
• Modesty on Limitations: The authors acknowledge that while the theoretical error estimates for non-convex domains suggest potential sub-optimality for $N > 2$ without specific regularity conditions, the numerical benchmarks with smooth solutions show that the method performs optimally in practice. They note that extensions to nonlinear equations, time-dependent problems, and other boundary conditions (Neumann, Robin) remain open challenges for future research. The work also suggests that the approach is adaptable to three-dimensional problems, citing related work in classical Lagrange finite elements.