Optimal boundary closures for diagonal-norm upwind SBP operators

By employing non-equispaced grid points near boundaries, boundary-optimized upwind finite-difference operators of orders up to nine are developed. The boundary closures are constructed within a diagonal-norm summation-by-parts (SBP) framework, ensuring linear stability on piecewise curvilinear multiblock grids. Boundary and interface conditions are imposed using either weak enforcement through simultaneous approximation terms (SAT) or strong enforcement via the projection method. The proposed operators yield significantly improved accuracy and computational efficiency compared with SBP operators constructed on equidistant grids. The resulting SBP–SAT and SBP–projection discretizations produce fully explicit systems of ordinary differential equations. The accuracy and stability properties of the proposed operators are demonstrated through numerical experiments for linear hyperbolic problems in one spatial dimension and for the compressible Euler equations in two spatial dimensions.

💡 Research Summary

This paper presents a systematic construction of diagonal‑norm upwind summation‑by‑parts (SBP) finite‑difference operators whose boundary closures are optimized for accuracy by employing non‑equispaced grid points near the domain boundaries. Traditional diagonal‑norm SBP operators, whether central or upwind, are usually built on uniform grids; consequently the formal order of accuracy drops near the boundaries (to at most p/2 for even p and (p‑1)/2 for odd p). This reduction limits the overall convergence rate and hampers long‑time simulations that demand high fidelity.

The authors address this limitation by allowing the first two grid points adjacent to each boundary to be placed at distances d₁h and d₂h (with h the interior spacing) rather than at uniform spacing. The parameters d₁ and d₂ become part of the design variables. Within the SBP framework, they define upwind operators D⁺ and D⁻ as
 D⁺ = H⁻¹(Q⁺ + B/2), D⁻ = H⁻¹(Q⁻ + B/2)
where H is a diagonal positive‑definite norm matrix, Q⁺ and Q⁻ are skew‑symmetric up to a negative‑semi‑definite matrix S = (Q⁺+Q⁺ᵀ)/2 = (Q⁻+Q⁻ᵀ)/2 ≤ 0, and B encodes the boundary contributions. The interior stencil of order p is obtained from a central (p+2)‑order narrow stencil by subtracting appropriately scaled forward/backward differences, yielding an off‑centered stencil that retains p‑th order accuracy while providing built‑in artificial dissipation through S.

The core of the construction is a constrained nonlinear optimization. The design variables include the entries of H (h₁,…,h_s), the entries of Q⁺ (the lower‑triangular part of the matrix), and the boundary spacings d₁, d₂. Constraints enforce: (i) the SBP symmetry conditions (Q⁺+Q⁺ᵀ=0, Q⁻+Q⁻ᵀ=0, Q⁺+Q⁻ᵀ=2S≤0); (ii) the required boundary accuracy, expressed via the error vectors e(q)=H x^q−1−(Q̂+B/2)x^q, which must vanish for q=1,…,p/2 (even p) or q=1,…,(p‑1)/2 (odd p). The objective function minimizes a linear combination of the leading‑order boundary error norms ‖e(p/2+1)‖₂ and ‖e(p/2+2)‖₂, thereby reducing the truncation‑error constant. The optimization is performed in Maple; d₁ and d₂ are restricted to