Nonsymmetric Interior Penalty Galerkin (NIPG)

• NIPG is a discontinuous Galerkin method defined by a nonsymmetric flux treatment that enhances computational efficiency and error control.
• It employs carefully chosen penalty parameters to ensure stability and optimal convergence rates in solving convection–diffusion, reaction–diffusion, and higher-order elliptic problems.
• Its implementation benefits from reduced element coupling costs and adaptive mesh strategies to effectively resolve complex geometries and boundary layer phenomena.

The Nonsymmetric Interior Penalty Galerkin (NIPG) method is a discontinuous Galerkin (DG) finite element technique for the numerical solution of partial differential equations (PDEs). Distinguished by its nonsymmetric flux treatment at element interfaces, NIPG is particularly relevant to convection–diffusion, reaction–diffusion, and higher-order elliptic problems, often on domains or problems exhibiting strong advection, boundary layers, or geometric complexity. NIPG is formulated by a specific choice of sign in flux terms, leading to an efficient, robust, and widely-studied DG approach, with distinct stability and error properties relative to symmetric (SIPG) or incomplete (IIPG) counterparts.

1. Mathematical Formulation

NIPG is characterized by the construction of a DG bilinear form with a nonsymmetric flux contribution. For an elliptic or parabolic PDE on a bounded polyhedral domain $\Omega\subset\mathbb{R}^d$ subdivided into a conforming mesh $\mathcal T_h$, the local broken finite element space is

$$V_h = \{v\in L^2(\Omega): v|_K\in \mathbb{P}_r(K),\ \forall K\in\mathcal{T}_h\}$$

with the set of (interior) faces $E_h$. On each face $e=K^+\cap K^-$, the DG formulation introduces the average and jump operators: $$\{u\} = \tfrac12(u^+ + u^-), \qquad [u] = u^- n^- + u^+ n^+,$$ where $n^\pm$ are the outward normals on $K^\pm$.

For scalar diffusion problems, the canonical NIPG bilinear form on $V_h$ is, for $\theta=-1$,

$$a_h(w,v) = \sum_{K}\int_{K}\nabla w \cdot \nabla v - \sum_{e}\int_{e}\{\nabla w\} \cdot [v] - \sum_{e}\int_{e} \{\nabla v\} \cdot [w] + \sum_{e}\int_{e} \frac{\sigma}{h_e} [w]\cdot [v].$$

For NIPG, the explicit minus sign in both flux terms yields a fully nonsymmetric operator, impacting coercivity and error structure (Kostas, 8 Nov 2025).

For multi-dimensional convection–diffusion or biharmonic problems, additional reaction, advection, and higher-order jump/penalty terms are included. The penalty parameter $\sigma$ is typically chosen as

$\sigma \ge C_0 r(r+1)$

with $C_0$ a geometric constant reflecting shape-regularity, and the penalty term is

$$\sum_{e\in E_h} \int_e \frac{\sigma}{h_e} [u_h]\cdot[v_h]\,ds.$$

For higher-order (e.g., biharmonic) problems on multipatch geometries, the NIPG bilinear form penalizes both function and derivative jumps, as detailed for isogeometric analysis (Moore, 2020).

2. Stability and Coercivity Properties

Stability of NIPG depends critically on penalty parameter magnitude and the underlying polynomial degree. Under standard mesh and solution regularity assumptions:

• For sufficiently large $\sigma$, coercivity is guaranteed:

$$a_h(v_h,v_h) \ge \beta \|v_h\|^2_{DG} - C \|v_h\|^2_{L^2(\Omega)}, \qquad \forall v_h\in V_h,$$

with $$\|v\|^2_{DG} = \sum_{K} \|\nabla v\|^2_{L^2(K)} + \sum_{e} (\sigma/h_e)\|[v]\|^2_{L^2(e)}$$.

• For time-dependent problems (e.g., reaction–diffusion or fractional PDEs), fully discrete NIPG schemes with backward Euler or L1-in-time and appropriate spatial penalties yield discrete energy inequalities and L^2-stability (Kostas, 8 Nov 2025, Maji et al., 2023).

The NIPG method typically requires a larger penalty ($\sigma$) to maintain coercivity compared with SIPG. For strong nonlinearities, penalty must be increased (by up to one order of magnitude) to preserve boundedness of discrete norms (Kostas, 8 Nov 2025).

On Shishkin/Bakhvalov layer-adapted meshes, penalty weights are varied according to region (smooth, boundary layer, internal interface) with sharply increased penalties in thin boundary/corner layers to preserve stability uniformly in $\varepsilon$ (Ma et al., 2023, Ma et al., 2023, Zhang et al., 2021).

For biharmonic isogeometric problems, coercivity is achieved for all nonnegative penalties in NIPG, with $\delta_0,\delta_1$ (penalty scaling constants) of order $p^2$ (polynomial degree), and boundedness in a mesh-dependent $|||v|||_h$ norm (Moore, 2020).

3. Error Estimates and Superconvergence

NIPG achieves optimal-order convergence in standard energy norms for piecewise polynomial degree $k$ or $r$:

• Elliptic and parabolic problems: In $L^2(\Omega)$, error is $O(h^{r+1})$, while in the DG-norm it is $O(h^r)$ (Kostas, 8 Nov 2025). For time-dependent, time-fractional problems, spatial error in the discrete energy-norm is $O(h^{k+1})$ (Maji et al., 2023).
• Singularly perturbed convection–diffusion: On Shishkin/Bakhvalov meshes, supercloseness of almost $k+1$ order is established between the NIPG solution and a layer-adapted interpolation (Radau+Lobatto in 1D, vertices-edge-element and $L^2$-projection in 2D) under the mesh-adapted energy norm. Postprocessing on macro-cells yields near-optimal $(k+1)$-th order global superconvergence in the energy-norm, confirmed numerically (Zhang et al., 2021, Ma et al., 2023, Ma et al., 2023).
• Higher-order (biharmonic) problems: In isogeometric contexts, NIPG achieves $O(h^{r-2})$ convergence rate in the mesh-dependent norm, with $r = \min(p+1,s)$ (B-spline degree, solution regularity) (Moore, 2020).

For time-fractional diffusion, superconvergent error $O(h^{k+1})$ is observed in a discrete energy-norm using Lobatto interpolants, while standard DG energy-norm only yields $O(h^k)$. Numerical experiments consistently confirm these results (Maji et al., 2023).

4. Implementation, Preconditioning, and Computational Aspects

NIPG offers efficiency in assembly and solution, especially exploiting its nonsymmetry:

• For standard DG, the nonsymmetric flux structure reduces the computational expense of element coupling compared to symmetric forms (Kostas, 8 Nov 2025).
• With reconstructed discontinuous approximation, arbitrary-order accuracy is achieved with one degree of freedom per element; this delivers significant savings in globally coupled unknowns for high-order accuracy. The reconstructed space is built by local least-squares fitting on patches, with the corresponding NIPG bilinear form imposed (Li et al., 2023).
• Preconditioning strategies exploit norm equivalence between the reconstructed high-order space and a piecewise-constant space, yielding a preconditioned system with mesh-independent condition number. Multigrid or algebraic multigrid can be used on the low-order preconditioner, resulting in Krylov iteration counts that are robust with respect to mesh and polynomial degree (Li et al., 2023).
• On graded or layer-adapted meshes (Shishkin, Bakhvalov), penalty weights are locally adapted and interpolation-based postprocessing is implemented to recover superconvergence while maintaining mesh-uniform stability (Zhang et al., 2021, Ma et al., 2023).
• In isogeometric multipatch settings, NIPG can be applied to NURBS-based H^2 spaces with non-matching grids, penalizing both function and gradient jumps to glue patches discretely (Moore, 2020).

5. Comparison with Symmetric and Incomplete IPDG Schemes

NIPG, SIPG, and IIPG are distinguished by the sign and completeness in their flux coupling terms ($\theta=1, -1, 0$):

• SIPG ($\theta=+1$): Fully symmetric, smallest required penalty for coercivity, most robust to strong nonlinearities. Yields bounded energy norms even for large coefficients. Slightly higher assembly cost (Kostas, 8 Nov 2025).
• NIPG ($\theta=-1$): Fully nonsymmetric, requires larger penalty for stability. If $\sigma$ is too small or nonlinearity is strong, energy norms may grow. Offers computational savings in assembly and is compatible with GMRES-type solvers (Kostas, 8 Nov 2025, Li et al., 2023).
• IIPG ($\theta=0$): Incompletely coupled, moderate stability and penalty requirements. Intermediate robustness (Kostas, 8 Nov 2025).

For high penalty parameter ($\sigma\gg1$), all three schemes converge to similar stability behavior. For moderate penalties, SIPG remains stable, IIPG retains partial robustness, while NIPG can destabilize unless penalties are increased (Kostas, 8 Nov 2025). This suggests that NIPG is best suited when computational throughput is prioritized or the nonsymmetry is advantageous (e.g., preconditioning). For strong nonlinear regimes or under-resolved grids, SIPG is generally preferred.

6. Practical Guidelines and Applications

The following practical recommendations have been established:

• Penalty parameter: Choose $\sigma\gtrsim C_0 r(r+1)$ (with $C_0\approx10$–$50$) and increase further for strong nonlinearities or advection-dominated problems. On nonuniform or layer-adapted meshes, increase penalty in regions corresponding to narrow layers or internal interfaces to maintain uniform stability (Kostas, 8 Nov 2025, Zhang et al., 2021, Ma et al., 2023).
• Polynomial degree: Select $r$ based on target spatial accuracy; NIPG matches SIPG in order but is more sensitive to under-penalization.
• Mesh design: For singularly perturbed problems, use Shishkin or Bakhvalov meshes to efficiently resolve internal layers. Adapt mesh transitions based on problem coefficients ($\varepsilon$, $\beta$) and use composite interpolation to preserve superconvergence (Zhang et al., 2021, Ma et al., 2023, Ma et al., 2023).
• Preconditioning: For large-scale high-order problems, construct low-order preconditioners using piecewise constant approximations; mesh-independence is achievable (Li et al., 2023).
• Verification: For each problem class, verify discrete $L^2$ norm boundedness via short, representative runs before large production computations (Kostas, 8 Nov 2025).
• Postprocessing: For optimal accuracy and superconvergence in boundary layers, implement macro-cell based higher-order interpolation postprocessing (Zhang et al., 2021, Ma et al., 2023).

NIPG is widely used in computational fluid dynamics, nonlinear optics, fractional diffusion, multipatch isogeometric analysis, and singular perturbation problems, often when interface flux upwinding or nonsymmetric coupling is advantageous. The method is robust and efficient for a wide problem class when its sensitivity to penalty parameters and stability considerations are properly managed.