🔢 mathematics

Penalty-scaling effects in nonsymmetric interior-penalty DG discretizations of viscous rotating shallow-water equations

This paper investigates the impact of interior-penalty parameter scaling on nonsymmetric interior-penalty DG discretizations of viscous rotating shallow-water equations, finding that while super-penalization ( $\beta=3$ ) can recover momentum $L^2$ accuracy, the standard scaling ( $\beta=1$ ) generally offers a superior accuracy-cost compromise for explicit implementations.

Original authors: Xue Zhang, Jingmin Xia, Xu Qian

Published 2026-07-20

📖 4 min read🧠 Deep dive

CC BY 4.0

Original authors: Xue Zhang, Jingmin Xia, Xu Qian

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the Earth's atmosphere and oceans as a giant, swirling dance floor. The dancers are huge masses of air and water, moving in complex patterns driven by the planet's spin and the sun's heat. To predict where a storm will go or how a current will shift, scientists use mathematical models called the "shallow-water equations." Think of these equations as the rulebook for the dance. However, real-world fluids aren't perfectly smooth; they have friction (viscosity) that slows things down, and they can get bumpy over underwater mountains.

To solve these rulebooks on a computer, scientists chop the dance floor into tiny puzzle pieces. This is where the "Discontinuous Galerkin" (DG) method comes in. Instead of forcing the dancers on neighboring puzzle pieces to hold hands perfectly, DG lets them step slightly out of sync, then uses a "penalty" to gently nudge them back into line. The big question is: how hard should you nudge them? If you nudge too lightly, the dancers might trip and the simulation crashes. If you nudge too hard, you might freeze the dance floor, making the computer calculation so slow it takes forever to finish. This paper investigates exactly how hard that nudge should be to keep the dance both accurate and fast.

The researchers, Xue Zhang, Jingmin Xia, and Xu Qian, tackled this problem using a specific version of the DG method called "Nonsymmetric Interior-Penalty Galerkin" (NIPG). They were studying the "viscous rotating shallow-water equations," which is a fancy way of saying they were modeling fluids that spin (like on Earth), have friction, and flow over a surface. The core of their study focused on a "penalty parameter," a dial that controls the strength of the nudge between puzzle pieces. They tested two main settings: a "standard" setting (where the nudge strength scales normally as the puzzle pieces get smaller) and a "super-penalized" setting (where the nudge gets incredibly strong, much stronger than usual).

Their findings reveal a fascinating trade-off, like choosing between a sports car and a tank. The "super-penalized" setting (where the penalty exponent $\beta = 3$ ) acts like a very strict referee. It does an excellent job of keeping the dancers (the momentum variables) perfectly in line, effectively eliminating any wobbly jumps between puzzle pieces. In fact, for the momentum part of the equation, this strict approach actually recovers the best possible mathematical accuracy. However, this strictness comes with a heavy price: it makes the computer's job incredibly stiff. The "stiffness" is like trying to run through waist-deep water; the computer has to take tiny, tiny steps in time to stay stable, which means the simulation runs much slower and costs significantly more computing power.

On the other hand, the "standard" setting (where the penalty exponent $\beta = 1$ ) is more like a relaxed coach. It allows for a bit more flexibility between the puzzle pieces. While it doesn't suppress the jumps as aggressively as the super-penalized version, the researchers found that it still provides excellent accuracy for the overall flow. Crucially, it doesn't slow down the computer nearly as much. In their tests, which included smooth rotating flows and waves moving over underwater mountains, the standard setting offered the best "bang for the buck." It gave them the right balance of accuracy and speed, whereas the super-penalized setting was overkill that made the calculations unnecessarily expensive without improving the final picture of the flow.

The paper also discovered something interesting about the different parts of the fluid model. While the super-penalized setting helped the "momentum" (the speed and direction of the flow), it didn't necessarily help the "geopotential" (which relates to the height of the water surface). In some cases, being too strict actually made the water height calculations slightly less accurate or no better at all. This suggests that for complex, real-world simulations involving rotating fluids, the "standard" approach is usually the smarter choice. It keeps the dance moving smoothly without freezing the dancers in place, proving that sometimes, a little bit of flexibility is better than a rigid, perfect alignment.

Technical Summary: Penalty-Scaling Effects in Nonsymmetric Interior-Penalty DG Discretizations of Viscous Rotating Shallow-Water Equations

Problem Statement
This study investigates the impact of penalty parameter scaling on the performance of Nonsymmetric Interior-Penalty Galerkin (NIPG) discretizations applied to the viscous rotating shallow-water equations (RSWEs) formulated in geopotential variables. While Discontinuous Galerkin (DG) methods are well-established for geophysical flows, the treatment of the elliptic viscous operator remains a critical design choice. Specifically, the paper addresses the trade-off between stability, accuracy, and computational stiffness introduced by the penalty law $\mu_e = \sigma h_e^{-\beta}$ . The authors compare the standard scaling ( $\beta = 1$ ) against a super-penalized scaling ( $\beta = 3$ ) to determine how these choices affect the coupled hyperbolic-parabolic system, particularly regarding momentum accuracy, geopotential behavior, and explicit time-step constraints.

Methodology
The authors develop a conservative DG formulation for the RSWEs where the hyperbolic terms are approximated using a local Lax–Friedrichs flux, and the viscous momentum diffusion is treated via a NIPG approach. The core of the methodology involves:

Analytical Framework: A linearized stability and error analysis is conducted around a smooth reference state. The authors establish consistency, continuity, and an exact coercivity identity for the diffusion form within a momentum DG seminorm. They derive a priori error estimates that distinguish between the momentum and geopotential variables, noting that the penalty exponent $\beta$ influences the convergence rates, particularly for even polynomial degrees.
Numerical Implementation: The semi-discrete system is advanced in time using a third-order Strong Stability Preserving Runge–Kutta (SSP-RK3) scheme. The implementation utilizes the Firedrake framework with the Unified Form Language (UFL).
Experimental Design: Four distinct numerical experiments are performed:
- A manufactured-solution convergence study to verify spatial rates.
- A penalty-sensitivity study varying $\beta \in \{1, 2, 3\}$ and the prefactor $\sigma$ to assess accuracy versus stiffness.
- A smooth rotating benchmark with non-zero Coriolis forcing to test geostrophic balance preservation.
- Topography-aware tests to verify well-balanced properties (lake-at-rest) and the evolution of small perturbations over non-flat bottoms.

Key Contributions
The paper makes four primary contributions:

Formulation: A conservative DG formulation in geopotential variables featuring a penalty-scaled NIPG momentum operator.
Theoretical Analysis: Explicit proofs of consistency, continuity, and a coercivity identity for the diffusion form in the momentum DG seminorm, alongside a conditional stability framework for the linearized coupled problem.
Error Estimation: Derivation of a priori error estimates that separate proven diffusion properties from assumptions on hyperbolic linearization, highlighting the conditional nature of optimal $L^2$ rates for NIPG.
Empirical Evidence: Comprehensive numerical evidence distinguishing the spatial accuracy, interface-jump control, and explicit computational stiffness associated with different penalty scalings.

Results
The numerical experiments yield several specific findings:

Convergence Rates: In manufactured-solution tests, super-penalization ( $\beta = 3$ ) successfully recovers the expected optimal $O(h^{k+1})$ momentum $L^2$ accuracy for odd polynomial degrees, whereas the standard scaling ( $\beta = 1$ ) often exhibits suboptimal rates for even degrees in the $L^2$ norm. However, the coupled geopotential variable does not necessarily exhibit the same improvement under super-penalization; in some cases, its accuracy degrades or remains unchanged.
Computational Stiffness: A significant trade-off is observed. Super-penalization ( $\beta = 3$ ) or large prefactors ( $\sigma$ ) induce severe algebraic stiffness, drastically reducing the admissible explicit time step (from $O(10^{-4})$ to $O(10^{-6})$ in tested cases) and increasing CPU cost.
Robustness and Balance: The scheme preserves the "lake-at-rest" equilibrium to machine precision regardless of penalty parameters. In rotating benchmarks, the scheme maintains geostrophic balance and captures gravity wave propagation without artificial damping, even under super-penalization.
Interface Control: While super-penalization significantly reduces interface jumps ( $J_{max}$ ), this mathematical control does not translate to improved overall solution accuracy in the tested regimes and comes at the cost of efficiency.

Significance and Claims
The paper claims that while super-penalization ( $\beta = 3$ ) offers theoretical benefits for restoring optimal elliptic $L^2$ approximation rates for momentum variables, it introduces prohibitive computational stiffness for explicit time integration in coupled geophysical systems. The authors conclude that for the explicit implementation considered, the standard scaling ( $\beta = 1$ ) generally provides a superior compromise between accuracy and computational cost. The study isolates the role of the penalty exponent, demonstrating that method parameters significantly determine the real operating range of a discretization. The authors suggest that future work should focus on semi-implicit time-integration methods to handle the stiffness introduced by NIPG terms and the extension to multi-layer systems for baroclinic dynamics.

More like this