High-order finite volume and discontinuous Galerkin methods are often stabilized by separate nonlinear devices for admissibility, entropy control, and oscillation suppression. This separation hides a simple geometric fact: all three act on the same cellwise candidate state. We propose a general framework (termed EPO) unifying fully discrete entropy stability, positivity/bound preservation, and spurious oscillation elimination. Starting from a candidate update, we scale along the ray anchored at its updated cell average. The admissible-state constraint, the entropy constraint, and the oscillation-suppressing constraint each define an admissibility radius on that ray, and the applied limiter is their minimum. The decisive analytical ingredient is a weak entropy stability at the level of the updated cell average. A two-point Lax–Friedrichs/Riemann-average entropy inequality yields local cell-average entropy budgets, and the same radial scaling mechanism behind Zhang–Shu positivity preservation lifts these weak budgets to strong quadrature-based entropy inequalities. The framework is therefore not a summation-by-parts, split-form, or flux-differencing construction: EPO acts on a candidate finite volume or discontinuous Galerkin update and converts weak average information into fully discrete nodal entropy stability. The construction also works for any prescribed finite family of convex entropy pairs. Each pair yields its own entropy radius, and taking the minimum enforces fully discrete entropy stability for all of them simultaneously. We prove the preservation of cell averages, invariant-set preservation, local and global strong entropy inequalities, stagewise budgets for strong-stability-preserving (SSP) Runge–Kutta methods, an SSP multistep variant that retains the designed high-order temporal accuracy, and extensions on rectangular and unstructured triangular meshes.

High-order accurate simulations of special relativistic hydrodynamics (RHD) are prone to numerical breakdown if intrinsic physical constraints (positive rest-mass density/pressure and subluminal velocity) are violated near strong discontinuities. In this work, we develop a robust and efficient physical-constraint-preserving (PCP) flux-limiting framework for high-order schemes, using finite-difference WENO as a representative example. By leveraging the geometric quasilinearization (GQL) representation, which equivalently reformulates the nonlinear RHD constraints into a family of linear inequalities, we integrate a Zalesak-type Flux-Corrected Transport (FCT) update into a scalar-style limiter that acts directly on conservative variables. A critical innovation is the explicit, non-iterative determination of limiting parameters via a rational stereographic parameterization of the GQL normal vector. This technique transforms the required worst-case minimization over auxiliary variables into a generalized Rayleigh-quotient formulation, allowing the optimal parameters to be obtained by solving small symmetric eigenvalue problems (2 × 2 in 1D; (d+1)×(d+1) in d dimensions). Relaxed variants are further introduced to reduce computational costs in multidimensions while retaining the PCP guarantee. Extensive numerical benchmarks ranging from 1D to 3D, including ultra-relativistic Riemann problems and astrophysical jets, demonstrate that the proposed method robustly enforces physical admissibility, sharply resolves discontinuities, and maintains design-order accuracy for smooth solutions.

Primitive variable recovery is a critical nonlinear subroutine in conservative relativistic hydrodynamics (RHD), often becoming a computational bottleneck for complex equations of state (EOS) in ultra-relativistic regimes. While analytical recovery formulas can be derived as closed-form algebraic identities, their floating-point evaluation may suffer from catastrophic cancellation when the Lorentz factor W >> 1; by contrast, generic Newton-type iterations may fail if the residual, admissible interval, and initial guess are not tied to a rigorous convergence argument. In this work, we first rewrite the primitive variable recovery problem for Synge-type EOSs in a normalized scalar form that separates the conservative-state kinematics from the EOS response, and introduce two canonical pressure residuals that are root-equivalent to the master equation. We then develop an EOS-shape-based endpoint-Newton framework with two complementary modes characterized by Type-L and Type-R EOS-shape conditions. The monotonicity, convexity, and one-sided initialization required by the Newton iterations are verified directly from the enthalpy law through shape functions. We show that the Taub-Mathews EOS satisfies the Type-L conditions and that the Ryu-Chattopadhyay-Choi EOS satisfies the Type-R conditions, each with equality in the corresponding curvature condition. This yields practical Newton solvers whose iterates converge monotonically to the physical pressure and achieve local quadratic convergence, without damping, ad hoc clipping, or case-specific fallback strategies. The numerical study combines baseline randomized benchmarks, theory-guided worst-case tests, and two-dimensional code-integration checks. The results show that the proposed recovery remains robust up to Lorentz factors W approximate to 10(6) for well-resolved states, while noticeable degradation occurs mainly near the floating-point admissibility frontier, where the conservative data themselves approach the limit of double-precision resolution.

This paper proposes an innovative, compact, characteristic-decomposition-free, and non-intrusive convex-oscillation-suppressing (COS) framework, termed COS(DG), for high-order discontinuous Galerkin (DG) discretizations of hyperbolic systems on general meshes. The method performs an entropy-guided convex scaling per time step, blending each DG polynomial with its cell average via coefficients determined by locally scale-and evolution-invariant entropy-induced distances. We prove that COS(DG) preserves the optimal convergence rate of the underlying DG scheme for smooth solutions, is L2-nonexpansive, and inherits entropy stability whenever the base DG formulation is entropy stable. The analysis bounds entropy-induced discrepancies of DG polynomials through L2 estimates and constructs a finite, compact, convex covering inside the admissible state set. We also provide a theoretical justification that the problem-independent COS procedure guarantees uniform oscillation suppression across diverse problems, thanks to the combined local scale-and evolution-invariant structures. The convex scaling strategy of COS(DG) offers a natural framework for integrating oscillation-suppressing and bound-preserving techniques. The COS approach is derivative-free, highly compact (using only immediate neighbors), and mode-independent (unified for modal and nodal DG), which facilitates parallelization and deployment on general meshes. A simple symmetry-preserving, scale-invariant entropy criterion detects interfaces near potential shocks and enables a local COS variant that further enhances resolution and reduces cost. Overall, COS(DG) offers nine salient advantages: (1) no characteristic decomposition; (2) provably optimal high-order convergence; (3) preservation of entropy stability (if present in the base DG method); (4) seamless integration with the Zhang-Shu bound-preserving limiter; (5) excellent compactness and easy parallelization using only neighboring data; (6) mode independence (modal/nodal unified); (7) simplicity via convex blending of solution and cell averages; (8) low overhead-applied only once per time step; and (9) local scale invariance (together with evolution invariance), en

We propose a bound-preserving (BP) Point-Average-Moment PolynomiAl-interpreted (PAMPA) scheme by blending third-order and first-order constructions. The originality of the present construction is that it does not need any explicit reconstruction within each element, and therefore the construction is very flexible. The scheme employs a classical blending approach between a first-order BP scheme and a high-order scheme that does not inherently preserve bounds. The proposed BP PAMPA scheme demonstrates effectiveness across a range of problems, from scalar cases to systems such as the Euler equations of gas dynamics. We derive optimal blending parameters for both scalar and system cases, with the latter based on the recent geometric quasi-linearization (GQL) framework of [Wu & Shu, SIAM Review, 65 (2023), pp. 1031–1073]. This yields explicit, optimal blending coefficients that ensure positivity and control spurious oscillations in both point values and cell averages. This framework incorporates a convex blending of fluxes and residuals from both high-order and first-order updates, facilitating a rigorous BP property analysis. Sufficient conditions for the BP property are established, ensuring robustness while preserving high-order accuracy. Numerical tests confirm the effectiveness of the BP PAMPA scheme on several challenging problems.

Constructing high-order schemes for ideal magnetohydrodynamics (MHD) that are simultaneously positivity-preserving (PP) and globally divergence-free (GDF) is challenging: conventional PP limiters typically break the interface continuity of the magnetic field required by GDF and, as shown in prior analysis (Wu, SIAM J. Numer. Anal., 2018 [1]), even small violations of the divergence constraint can destroy positivity. A straightforward workaround to make PP limiting compatible with GDF is to add a corrective amount to the total energy; however, this sacrifices energy conservation. This paper proposes a discontinuous Galerkin (DG) approach that couples a novel convex oscillation-suppressing (COS) procedure with two conditionally conservative energy-modification strategies (an optimization-based limiter and a closed-form energy-scaling limiter) to enforce the pointwise positivity without altering the GDF magnetic field. We show that energy-only conservative PP enforcement is feasible if and only if a computable condition holds; otherwise, we apply a local energy-compensation fallback. This condition fails only in a negligible fraction of cells: typically below 0.03% and never exceeding 1% in our tested cases. The COS mechanism damps spurious oscillations via a cell-wise entropy-Hessian-weighted extension operator, avoiding derivative jumps and additional tunable parameters beyond a single problem-independent factor. To establish a rigorous PP property, we introduce auxiliary magnetic-field cell averages updated in the same manner as other conservative variables and employ an optimal convex decomposition on Gauss-Lobatto boundary nodes, which also relaxes CFL constraints and reduces the number of PP correction points. Extensive tests-including the Orszag-Tang vortex, rotor, strong blasts with ambient plasma-beta of about 2.5 & times; 10-6, shock-cloud interaction, and a Mach-10000 jet-demonstrate robustness, high-order accuracy, preservation of the GDF property, and near-conservative behavior of the global total energy delivered by the proposed optimization and scaling energy-modification techniques.

We derive an explicit discrete energy identity for rational time discretizations generated by the first-subdiagonal Padé approximants of the exponential for solving linear seminegative problems. This work extends the diagonal Padé energy laws in [Z. Sun, Y. Wei, and K. Wu, SIAM J. Numer. Anal., 60 (2022)] to the first-subdiagonal family. The main new ingredient is an explicit Cholesky-type factorization of the energy coefficient matrix associated with the semi-inner-product terms in the discrete energy identity. The construction and proof of this factorization are nontrivial, since the matrix entries are alternating sums of Padé coefficients and the triangular factor has a parity-dependent factorial structure. We prove the factorization by reducing it to scalar rational identities and establishing them through finite product reductions and telescoping summations. Together with a β-coefficient cancellation, the factorization yields an exact discrete energy law that recovers the classical unconditional contractivity for linear seminegative problems. Numerical experiments adapted from the diagonal Padé energy-law setting illustrate the predicted order and verify the discrete dissipation identity.

A recurring challenge in science and engineering is the model-reality gap, where trusted legacy simulators lose fidelity due to unresolved physics or structural incompleteness. This challenge has motivated remedies ranging from imperfect mechanistic models to fully data-driven surrogates. Here, we address this gap with Alternating Neural Integrators (ANI), a non-intrusive reuse-and-correct framework for upgrading executable legacy simulators without requiring access to or modification of their internal implementation. Guided by operator-splitting principles, ANI alternates the evolution of a fixed prior simulator with a learned neural correction that targets structured discrepancy between the prior and the supervisory data. We show that ANI can recover missing coupling in chaotic systems and act as an effective subgrid correction in turbulence, improving dynamical fidelity where prior models drift. In selected mechanistically structured settings, post hoc symbolic distillation yields compact hypotheses and, in controlled benchmarks, supports closed-loop refinement of the prior. By combining data-driven flexibility with reusable scientific simulators in a theoretically grounded framework, this work provides a practical gray-box route for systematically upgrading existing computational infrastructure when a callable prior and supervisory data are available.

This paper establishes the first rigorous superconvergence theory for semidiscrete and fully discrete central discontinuous Galerkin (CDG) methods for linear hyperbolic equations on overlapping meshes. While the optimal $L^2$ convergence of $\mathbb{Q}^k$ CDG schemes was established on uniform Cartesian meshes by Liu, Shu, and Zhang [ SIAM J. Numer. Anal.}, 56 (2018), pp. 520--541], their observed $\mathcal{O}(h^{k+2})$ pointwise superconvergence has remained unproven, due to the loss of standard single-mesh Galerkin orthogonality inherent in the CDG overlapping structure. To overcome this fundamental barrier, we introduce a projection-correction framework that identifies a hidden superconvergent mechanism: an asymptotic weak residual cancellation in one dimension, and a high-order cancellation-by-aggregation (HOCA) mechanism in multiple dimensions. This HOCA approach overcomes the analytical challenge posed by coupled primal-dual directional residuals, recovering critical error cancellation properties absent from the standard variational formulation. Consequently, we provide the rigorous proof of the conjectured $\mathcal{O}(h^{k+2})$ pointwise superconvergence in the discrete $\ell^{\infty}$ norm across all superconvergent points. Furthermore, we reveal that under a systematically corrected initialization, this framework yields a previously undiscovered, stronger cell-average superconvergence estimate of order $\mathcal{O}(h^{\min\{2k+1,k+3\}})$. The theory is extended to fully discrete explicit Runge--Kutta CDG schemes, where stagewise corrected errors are constructed to preserve spatial superconvergence up to temporal truncation errors, yielding a stable reconstruction-based postprocessing estimate. Numerical experiments in one and two spatial dimensions confirm the sharpness of the theoretical rates.

This paper proposes and analyzes the OECDG method, a novel central discontinuous Galerkin (CDG) scheme that integrates the strengths of the oscillation-eliminating (OE) approach [M. Peng, Z. Sun, and K. Wu, Math. Comp., 94 (2025), pp. 1147--1198] within the CDG framework for general hyperbolic conservation laws. The OECDG method incorporates a new OE procedure, along with an innovative dual damping mechanism, to enhance oscillation control in CDG schemes. Unlike OE procedures in standard DG methods that use intercell jumps as a modal filter, our new OE mechanism draws inspiration from CDG-based numerical dissipation and employs a convex blending strategy, leveraging overlapping solutions to enhance stencil compactness and resolution without characteristic decomposition. We rigorously derive optimal error estimates for the fully discrete OECDG method through several key theoretical advancements, filling a gap in the error analysis of fully discrete CDG schemes, including original linear CDG methods without oscillation control. First, we establish the approximate skew-symmetry and weak boundedness of the CDG spatial discretization, which underpins the fully discrete stability analysis. Utilizing these properties, we then prove the linear stability of CDG methods coupled with Runge--Kutta time discretization via matrix transfer techniques. These foundational results enable us to derive fully discrete optimal error estimates for the OECDG method---a challenging task due to the method's nonlinear nature, even for linear advection equations. Extensive numerical experiments validate the theoretical findings and demonstrate the efficacy of the OECDG method across a variety of hyperbolic conservation laws, including linear and nonlinear problems such as the convection equation, the Burgers equation, a traffic flow model, and Euler equations. The results confirm the OECDG method's ability to achieve optimal convergence rates, robustly eliminate spurious oscillations, and accurately capture complex wave structures across diverse test cases.

This paper proposes high-order accurate bound-preserving (BP) finite volume methods on adaptive moving structured meshes for two- and three-dimensional special relativistic hydrodynamics (RHD). The BP property here includes the positivity of rest-mass density and pressure, the subluminal constraint on fluid velocity, as well as the minimum entropy principle established in [1]. The methods are built on the time-dependent coordinate transformation from the computational domain to the physical domain, appropriate discretization of the geometric conservation laws (GCLs), a global Lax-Friedrichs (LF) type numerical flux incorporating the mesh metrics, and the explicit strong-stability-preserving Runge-Kutta time discretizations. Preserving the minimum entropy principle is nontrivial, as the commonly used LF splitting property no longer holds in general. To address this, a weak LF splitting property, compatible with the minimum entropy principle, is introduced. A rigorous BP analysis is conducted based on the weak LF splitting property, the discrete GCLs, and the geometric quasilinearization (GQL) approach in [2, 3]. Finally, various numerical examples in two and three dimensions are presented to validate the high-order accuracy, high resolution, efficiency, and BP property of the proposed methods.

The PAMPA (Point-Average-Moment PolynomiAl-interpreted) method, proposed in [1], is a compact active-flux-type framework that combines conservative and nonconservative formulations of hyperbolic problems. In this paper, we develop a novel positivity-preserving (PP) PAMPA scheme for the ideal magnetohydrodynamics (MHD) equations on Cartesian grids. The method enforces a point-value-level discrete divergence-free (DDF) constraint on the intermediate interface states used in the update procedure, while the fully evolved continuous quadratic representation is not claimed to be exactly divergence-free at every instant. Building on our recent one-dimensional invariant-domain-preserving PAMPA framework [2], we extend the methodology to the multidimensional MHD system. The proposed scheme features an automatic PP update of interface point values through a new nonconservative reformulation, together with a local DDF projection. The cell-average update is provably PP under a mild a priori positivity condition on a single cell-centered value and combines four ingredients: (i) a DDF constraint at interface point values, (ii) a PP limiter applied only to the cell-centered value, (iii) a PP numerical flux with properly estimated wave speeds, and (iv) a suitable discretization of the Godunov–Powell source term. The PP proof for cell averages is carried out within the geometric quasi-linearization (GQL) framework [3], which transforms the nonlinear pressure-positivity constraint into an equivalent linear form. The resulting scheme avoids explicit polynomial reconstruction in practice, is compatible with arbitrarily high-order strong-stability-preserving time discretizations, and is straightforward to implement. To enhance robustness and resolution, we introduce a problem-independent troubled-cell indicator based on a Lax-type entropy criterion that uses only two characteristic speeds computed from cell averages and preserves symmetry, together with a convex oscillation suppression (COS) mechanism with a new distance measuring solution differences between target and neighboring cells. Extensive tests, including a rotated shock tube, a blast wave with plasma beta as low as 2.51×10−6, and a jet with Mach number up to 104, demonstrate high-order accuracy, sharp resolution of complex MHD structures, strong-shock robustness, and bounded analytical/discrete divergence diagnostics when these are evaluated from the unlimited representation. To the best of our knowledge, this is the first active-flux-type method for ideal MHD that combines a rigorous PP proof for cell averages, automatic admissibility of interface point values, and a point-value-level DDF mechanism in the update procedure.

A discrete entropy inequality is the principal nonlinear stability estimate available for systems of conservation laws, and evaluating it presupposes a physically admissible state. So far, however, the two have been secured separately. Entropy-stable schemes are almost always semi-discrete, are built around one selected entropy pair, and take for granted the positivity of density and pressure that makes the entropy well defined in the first place, whereas bound-preserving limiters keep the solution admissible but deliver no entropy estimate. For the special relativistic Euler equations, the two cannot be separated at all, since the conservative-to-primitive map is implicit, and an inadmissible state therefore has no entropy to correct. Here we construct high-order discontinuous Galerkin and finite volume schemes that, to our knowledge, for the first time, are entropy stable in the fully discrete sense for an arbitrary prescribed finite family of convex entropy pairs, a property we call multi-entropy stability, and are provably admissible wherever an entropy is evaluated. All of this is achieved by a single cellwise projection, and neither conservation nor the design order is lost. The construction rests on relativistic causality, which bounds every characteristic speed by the speed of light. Consequently, the numerical viscosity can be fixed once for all states and all equations of state, and one two-point building block then serves the whole entropy family. Since only the convexity of the admissible set and this speed bound are used, the same route remains open for related systems. Finally, in computations with four equations of state, the schemes retain high-order accuracy close to vacuum, produce no inadmissible state in strong shocks, near-vacuum shock–vortex interaction or jets with Lorentz factor above 70, and confirm the monotone decay of every enforced discrete entropy.

This paper explores numerical schemes for Temple-class systems, which are integral to various applications including one-dimensional two-phase flow, elasticity, traffic flow, and sedimentation. Temple-class systems are characterized by conservative equations, with different pressure function expressions leading to specific models such as the Aw-Rascle-Zhang (ARZ) traffic model and the sedimentation model. Our work extends existing studies by introducing a moving mesh approach to address the challenges of preserving non-convex invariant domains, a common issue in the numerical simulation of such systems. Our study outlines a novel bound-preserving (BP) and conservative numerical scheme, designed specifically for non-convex sets in Temple-class systems, which is critical for avoiding non-physical solutions and ensuring robustness in simulations. We develop both local and global BP methods based on finite difference schemes, with numerical experiments demonstrating the effectiveness and reliability of our methods. Furthermore, a parameterized flux limiter is introduced to restrict high-order fluxes and maintain bound preservation. This innovation marks the first time such a parameterized approach has been applied to non-convex sets, offering significant improvements over traditional methods. The findings presented extend beyond theoretical implications, as they are applicable to general Temple-class systems and can be tailored to ARZ traffic flow networks, highlighting the versatility and broad applicability of our approach. The paper contributes significantly to the field by providing a comprehensive method that maintains the physical and mathematical constrains of Temple-class systems.

Admissible states in hyperbolic systems and related equations often form a convex invariant domain. Numerical violations of this domain can lead to loss of hyperbolicity, resulting in illposedness and severe numerical instabilities. It is therefore crucial for numerical schemes to preserve the invariant domain to ensure both physically meaningful solutions and robust computations. For complex systems, constructing invariant-domain-preserving (IDP) schemes is highly nontrivial and particularly challenging for high-order accurate methods. This paper presents a comprehensive survey of IDP schemes for hyperbolic and related systems, with a focus on the most popular approaches for constructing provable IDP schemes. We first give a systematic review of the fundamental approaches for establishing the IDP property in first-order accurate schemes, covering finite difference, finite volume, finite element, and residual distribution methods. Then we focus on two widely used and actively developed classes of high order IDP schemes as well as their recent developments, most of which have emerged in the past decade. The first class of methods seeks an intrinsic weak IDP property in high-order schemes and then designs polynomial limiters to enforce a strong IDP property at the points of interest. This generic approach applies to high-order finite volume and discontinuousGalerkin schemes. The second class is based on the flux limiting approaches, which originated from the flux-corrected transport method and can be adapted to a broader range of spatial discretizations, including finite difference and continuous finite element methods. In this survey, we elucidate the main ideas in the construction of IDP schemes, provide some new perspectives and insights, with extensive examples, and numerical experiments in gas dynamics and magnetohydrodynamics.

One major challenge in developing accurate and robust numerical schemes for compressible Euler equations arises due to the potential emergence of discontinuous structures in the solution. Recently proposed low-dissipation central-upwind (LDCU) schemes achieve sharp resolution of these structures without introducing spurious oscillations. However, unlike many other Godunov-type methods, the LDCU schemes cannot be written as a convex combination of first-order positivity-preserving (PP) schemes. Therefore, the PP property of the LDCU schemes cannot be analyzed by standard techniques. In this paper, we overcome this difficulty by first decomposing the studied schemes into a convex combination of several intermediate solution states, and then analyzing their PP properties. The performed analysis helps us to construct PPLDCU schemes for Euler equations of compressible gas dynamics, guaranteeing the positivity of computed density and pressure. To achieve the PP property, the built-in anti-diffusion terms in the two-dimensional case and the piecewise linear reconstruction procedure in both the one- and two-dimensional cases are redesigned. The effectiveness and robustness of the proposed PPLDCU schemes are demonstrated in several challenging numerical examples.

This paper proposes and analyzes a class of essentially non-oscillatory central discontinuous Galerkin (CDG) methods for general hyperbolic conservation laws. First, we introduce a novel compact, non-oscillatory stabilization mechanism that effectively suppresses spurious oscillations while preserving the high-order accuracy of CDG methods. Unlike existing limiter-based approaches that rely on large stencils or problem-specific parameters for oscillation control, our dual damping mechanism is inspired by CDG-based numerical dissipation and leverages overlapping solutions within the CDG framework, significantly enhancing stability while maintaining compactness. Our approach is free of problem-dependent parameters and complex characteristic decomposition, making it efficient and robust. Second, we provide a rigorous stability and optimal error analysis for fully discrete Runge-Kutta (RK) CDG schemes, addressing a gap in the theoretical understanding of these methods. Specifically, we establish the approximate skew-symmetry and weak boundedness of the CDG discretization. These results enable us to rigorously analyze the fully discrete error estimates for our oscillation-eliminating CDG (OECDG) method, a challenging task due to its nonlinear nature, even for linear advection equations. Building on this framework, we reformulate nonlinear oscillation-eliminating CDG schemes as linear RK CDG schemes with a nonlinear source term, extending error estimates beyond the linear case to schemes with nonlinear oscillation control. While existing error analyses for DG or CDG schemes have largely been restricted to linear cases without nonlinear oscillation-control techniques, our analysis represents an important theoretical advancement. Experiments validate the theoretical findings and demonstrate the effectiveness of the OECDG method.

Quadrature-based moment methods (QBMM) provide tractable closures for multiscale kinetic equations, with diverse applications across aerosols, sprays, and particulate flows, etc. However, for the derived hyperbolic moment-closure systems, seeking numerical schemes preserving moment realizability is essential yet challenging due to strong nonlinear coupling and the lack of explicit conservative-to-flux maps. This paper proposes and analyzes a provably realizability-preserving finite-volume method for five-moment systems closed by the two-node Gaussian-EQMOM and three-point HyQMOM. Rather than relying on kinetic fluxes, we recast the realizability condition into a nonnegative quadratic form in the moment vector, reducing the original nonlinear constraints to bilinear inequalities amenable to analysis. On this basis, we construct a tailored Harten–Lax–van Leer (HLL) flux with rigorously derived wave speeds and intermediate states that embed realizability directly into the flux evaluation. We prove sufficient realizability-preserving conditions under explicit Courant–Friedrichs–Lewy (CFL) constraints in the collisionless case, and for BGK relaxation, we obtain coupled time-step conditions involving a realizability radius; a semi-implicit BGK variant inherits the collisionless CFL. From a multiscale perspective, the analysis yields stability conditions uniform in the relaxation time and supports stiff-to-kinetic transitions. A practical limiter enforces strict realizability of reconstructed interface states without degrading accuracy. Numerical experiments demonstrate the accuracy, robustness in low-density regions, and realizability for both closures. This framework unifies realizability preservation for solving hyperbolic moment systems with complex closures and extends naturally to higher-order space–time discretizations.

This paper develops high-order accurate, well-balanced (WB), and positivity-preserving (PP) finite volume schemes for shallow water equations on adaptive moving structured meshes. The mesh movement poses new challenges in maintaining the WB property, which not only depends on the balance between flux gradients and source terms but is also affected by the mesh movement. To address these complexities, the WB property in curvilinear coordinates is decomposed into flux source balance and mesh movement balance. The flux source balance is achieved by suitable decomposition of the source terms, the numerical fluxes based on hydrostatic reconstruction, and appropriate discretization of the geometric conservation laws (GCLs). Concurrently, the mesh movement balance is maintained by integrating additional schemes to update the bottom topography during mesh adjustments. The proposed schemes are rigorously proven to maintain the WB property by using the discrete GCLs and these two balances. We provide rigorous analyses of the PP property under a sufficient condition enforced by a PP limiter. Due to the involvement of mesh metrics and movement, the analyses are nontrivial, while some standard techniques, such as splitting high-order schemes into convex combinations of formally first-order PP schemes, are not directly applicable. Various numerical examples validate the high-order accuracy, high efficiency, WB, and PP properties of the proposed schemes.

This paper establishes the minimum entropy principle (MEP) for the relativistic Euler equations with a broad class of equations of state (EOSs) and addresses the challenge of preserving the local version of the discovered MEP in high-order numerical schemes. At the continuous level, we find out a family of entropy pairs for the relativistic Euler equations and provide rigorous analysis to prove the strict convexity of entropy under a necessary and sufficient condition. At the numerical level, we develop a rigorous framework for designing provably entropy-preserving high-order schemes that ensure both physical admissibility and the discovered MEP. The relativistic effects, coupled with the abstract and general EOS formulation, introduce significant challenges not encountered in the nonrelativistic case or with the ideal EOS. In particular, entropy is a highly nonlinear and implicit function of the conservative variables, making it particularly difficult to enforce entropy preservation. To address these challenges, we establish a series of auxiliary theories via highly technical inequalities. Another key innovation is the use of geometric quasi-linearization (GQL), which reformulates the nonlinear constraints into equivalent linear ones by introducing additional free parameters. These advancements form the foundation of our entropy-preserving analysis. We propose novel, robust, locally entropy-preserving high-order frameworks. A central challenge is accurately estimating the local minimum of entropy, particularly in the presence of shock waves at unknown locations. To address this, we introduce two new approaches for estimating local lower bounds of specific entropy, which prove effective for both smooth and discontinuous problems. Numerical experiments demonstrate that our entropy-preserving methods maintain high-order accuracy while effectively suppressing spurious oscillations.