New Adaptive Numerical Methods Based on Dual Formulation of Hyperbolic Conservation Laws

Abstract: In this paper, we propose an adaptive high-order method for hyperbolic systems of conservation laws. The proposed method is based on a dual formulation approach: Two numerical solutions, corresponding to conservative and nonconservative formulations of the same system, are evolved simultaneously. Since nonconservative schemes are known to produce nonphysical weak solutions near discontinuities, we exploit the difference between these two solutions to construct a smoothness indicator (SI). In smooth regions, the difference between the conservative and nonconservative solutions is of the same order as the truncation error of the underlying discretization, whereas in nonsmooth regions, it is ${\cal O}(1)$. We apply this idea to the Euler equations of gas dynamics and define the SI using differences in the momentum and pressure variables. This choice allows us to further distinguish neighborhoods of contact discontinuities from other nonsmooth parts of the computed solution. The resulting classification is used to adaptively select numerical discretizations. In the vicinities of contact discontinuities, we employ the low-dissipation central-upwind numerical flux and a second-order piecewise linear reconstruction with the slopes computed using an overcompressive SBM limiter. Elsewhere, we use an alternative weighted essentially non-oscillatory (A-WENO) framework with the central-upwind finite-volume numerical fluxes and either unlimited (in smooth regions) or Ai-WENO-Z (in the nonsmooth regions away from contact discontinuities) fifth-order interpolation. Numerical results for the one- and two-dimensional compressible Euler equations show that the proposed adaptive method improves both the computational efficiency and resolution of complex flow features compared with the non-adaptive fifth-order A-WENO scheme.

• The paper introduces a dual formulation that solves both conservative and primitive representations concurrently to establish a robust smoothness indicator.
• It employs adaptive high-order schemes that switch between fifth-order and low-dissipation solvers based on local smoothness criteria.
• It demonstrates sharper resolution of flow features and improved computational efficiency compared to standard A-WENO methods in various benchmark tests.

Adaptive High-Order Schemes via Dual Formulation for Hyperbolic Conservation Laws

Introduction and Motivation

The paper "New Adaptive Numerical Methods Based on Dual Formulation of Hyperbolic Conservation Laws" (2601.20000) develops a novel adaptive high-order methodology for the numerical solution of hyperbolic systems of conservation laws. The approach is structurally centered on a dual formulation (DF), where both conservative and nonconservative (primitive-variable) representations of the governing equations are solved concurrently. This dual evolution furnishes a robust adaptive criterion derived from the local discrepancy between the two solutions, capitalizing on the recognized deficiency of nonconservative schemes in the vicinity of discontinuities to design a smoothness indicator (SI). Consequently, the scheme can locally select among high-order and low-dissipation solvers in a region-aware fashion, optimizing both resolution and computational efficiency.

Dual Formulation and Smoothness Indicator Construction

Let $U_t + F(U)_x + G(U)_y = 0$ define a generic $d$-dimensional hyperbolic conservation law, with conserved variables $U$. The associated nonconservative (primitive variable) form is

$$V_t + \widetilde{F}(V)_x + \widetilde{G}(V)_y = B(V) V_x + C(V) V_y,$$

where $V$ denotes primitive variables, and $B(V),\,C(V)$ encode the nonconservative products. For smooth solutions, both forms are equivalent, but near discontinuities, nonconservative solvers are well known to yield nonphysical results due to the incorrect weak solution selection.

In the proposed DF scheme, both formulations are solved at each time step. The key innovation is defining the SI through the cellwise difference between the conservative and nonconservative solutions, projected onto selected physically relevant variables (e.g., momentum, pressure for compressible Euler equations). Specifically, in smooth regions, this difference scales with the truncation error, while in nonsmooth regions, it is $\mathcal{O}(1)$. The SI is then locally analyzed to classify cells as:

• S (Smooth): Difference below threshold; high-order interpolation is admissible.
• RC (Contact): Momentum difference large, pressure difference small; indicates contact discontinuities.
• RNC (Rough Non-Contact): Both momentum and pressure differences large; corresponds to shocks and other discontinuities.

Thresholds ($\kappa_{\rho u}, \kappa_p, \ldots$) are tuned per test, but the method's robustness with respect to threshold and mesh refinement is demonstrated.

Adaptive Scheme Design

The fully discrete algorithm advances the split solution using different discretizations, depending on SI-based local classification:

• Region S: Fifth-order A-WENO with central-upwind (CU) flux and unlimited interpolation.
• Region RNC: Fifth-order A-WENO, but with Ai-WENO-Z nonlinear interpolation in the variable's characteristic fields.
• Region RC: Second-order low-dissipation central-upwind (LDCU) scheme with an overcompressive (SBM) slope limiter, also in the characteristic fields.

The primitive (nonconservative) system is evolved with a simplified fifth-order A-WENO discretization, used only to compute the SI; its result is overwritten each time step by the conservative solution. Notably, for the Euler equations, the SI employs the difference between conservative and nonconservative momenta and pressures, leveraging the continuity of pressure across contacts.

All schemes employ third-order SSP Runge-Kutta for time integration, with a typical CFL number of $0.45$.

Numerical Results

One-Dimensional Benchmarks

Shock-Density Wave Interaction: The adaptive scheme matches or slightly outperforms A-WENO in accuracy on the same mesh, with significant CPU time savings. For equal computational budgets, the adaptive method provides superior resolution of fine structures.

Shock-Entropy Wave Interaction: The adaptive method provides sharper results for high-frequency entropy waves compared to standard A-WENO at equal computational cost.

Blast Wave: The method yields significantly sharper resolution of contact waves, outperforming A-WENO, due to correct RC region detection and activation of low-dissipation/slope-overcompressive methods near contacts.

Two-Dimensional Benchmarks

2D Riemann Problems: For both canonical test configurations, the adaptive approach yields more accurate shock and vortex structures and captures Kelvin-Helmholtz roll-up with less numerical dissipation than A-WENO. The region classifier reliably tracks contact and shock neighborhoods, even under mesh refinement and threshold variation.

Implosion Problem: The adaptive method more sharply resolves jet-like structures and high-gradient interfaces, outperforming A-WENO on identical meshes.

Rayleigh-Taylor Instability: The adaptive method captures significantly finer-scale structures and Kelvin-Helmholtz vortices in the mixing region, which are smeared or missing in the non-adaptive solvers.

In all cases, the computational overhead associated with region classification and dual solves is more than compensated by savings from restricting costly WENO reconstruction to rough regions only.

Theoretical and Practical Implications

The proposed DF-based adaptivity successfully leverages the inherent inconsistency of nonconservative discretizations near discontinuities as a signal for adaptivity, rather than a pathology to be suppressed. This results in an SI that is more physically discriminative, allowing for region-specific numerical treatment, which is particularly beneficial for complex multi-structure flows where both accuracy and efficiency are critical.

Practical implications include:

• Enhanced Resolution: Selective application of low-dissipation/high-sharpness methods near contacts and high-order methods elsewhere produces both sharper and more efficient results.
• Computational Efficiency: The method achieves higher-fidelity results at the same or reduced computational cost.
• Robust Region Classification: Demonstrated stability to variation in SI thresholds and mesh refinement.
• Framework Flexibility: The DF paradigm allows for modular incorporation of advanced high-order, adaptive, and structure-preserving schemes.

Theoretical implications extend to error analysis of hybrid algorithms, pointing toward new frameworks for SI design in complex, multi-regime systems, including multifluid and multiphase flows, as well as for asymptotic-preserving schemes.

Future Directions

The authors suggest potential extensions to compressible multifluid models, as the DF approach is naturally compatible with systems formulated in both conservative and primitive variables (e.g., five- or seven-equation models). Integration with well-balanced, positivity-preserving, and bound-preserving techniques is natural. Furthermore, the approach's modular SI design promises direct applicability to locally adaptive mesh refinement and anisotropic adaptivity.

Conclusion

The dual formulation-based adaptive numerical framework introduced in this work delivers an effective synthesis of region-aware adaptivity and high-order accuracy for hyperbolic conservation laws (2601.20000). By utilizing the discrepancy between conservative and nonconservative solutions as a robust smoothness indicator, the method systematically outperforms standard non-adaptive high-order schemes, both in efficiency and in the resolution of flow features such as contact discontinuities and shock-dominated structures. This DF-based paradigm represents a concrete step toward general, efficient, and reliable simulation strategies for complex compressible flows involving intricate interactions between smooth and nonsmooth structures. Future research will likely address the algorithmic extension to multifluid and asymptotic-preserving contexts and deeper theoretical analysis of SI-driven adaptivity.