Preserving the admissible set of the ideal magnetohydrodynamics (MHD) equations is important not only for producing physically meaningful numerical solutions, but more importantly for achieving robust computations. In this paper, we develop an optimization-based limiter to enforce admissibility while preserving global conservation and accuracy. For an easy and efficient projection, we decompose the admissible set into slices parameterized by the magnetic energy, so that the MHD projection reduces to a one-dimensional minimization, which can be solved efficiently by the Brent method. The splitting method can be used to efficiently solve the global minimization problem of the optimization-based limiter, which can be used to enforce cell average admissibility in discontinuous Galerkin (DG) schemes, and pointwise admissibility can be further enforced by the Zhang-Shu positivity-preserving limiter. We apply the limiter to high-order DG schemes and present numerical results for a few representative MHD problems.

Finite difference weighted essentially non-oscillatory (WENO) schemes fill an important need in computational science—in multi-dimensions they can carry out high accuracy simulations at a fraction of the cost of finite volume and discontinuous Galerkin alternatives. Alternative finite difference WENO (AFD-WENO) schemes are variants of WENO schemes which provide the added flexibility of handling structured meshes with complex geometry. AFD-WENO schemes have also been designed recently to handle hyperbolic PDEs with non-conservative products (Balsara et al. in Commun Appl Math Comput, 2024. 10.1007/s42967-024-00374-1), thereby dramatically increasing the range of application areas to which they can be applied. However, one of the deficiencies of higher order WENO schemes stems from the fact that there are stringent problems where the schemes are not Physical Condition Preserving (PCP). This paper rectifies the above-mentioned deficiency by designing and documenting AFD-WENO schemes with a PCP property. The essential idea is to selectively hybridize a low order scheme, which is PCP, with a high order scheme, which may sometimes lose the PCP property. These ideas are packaged into a formulation that works with the fluctuation form for hyperbolic systems, and yet, our formulation is as easy to implement as PCP conditions for conservation laws. The formulation can work seamlessly for conservation laws as well as PDEs with non-conservative products. When the PDE has a conservation form in some limits, the formulation is designed to retrieve that conservation form. We numerically observe that the formulation does not damage the high order accuracy for problems with smooth flow. Several extremely stringent test problems are presented in this paper to illustrate the value of the PCP methods developed here.

Weighted essentially nonoscillatory (WENO) schemes are a popular class of numerical methods for solving hyperbolic conservation laws. Since WENO schemes are designed to deal with problems with both complicated solution structures and discontinuities / sharp gradient regions, their sophisticated nonlinear properties and high-order accuracy require more operations than many other schemes. The methodology of hybrid methods is an effective approach to decrease the computational costs and dissipation errors of WENO schemes and achieve better resolution. One of the key components for the success of hybrid WENO schemes is the application of a robust and efficient troubled-cell indicator, which detects the computational cells where the solution loses regularity. Recently, troubled-cell indicators based on artificial neural networks (ANNs) have been developed in the literature, which have the advantage of less dependence on tunable parameters and being more robust than many traditional troubled-cell indicators, and such ANN based troubled-cell indicators have been applied to hybrid finite difference WENO schemes effectively. Motivated by these works, in this paper we develop a hybrid finite volume WENO method with an ANN based troubled-cell indicator for solving hyperbolic conservation laws. While the finite difference WENO schemes are more efficient than the finite volume WENO schemes for multidimensional problems on uniform grids, the finite volume WENO schemes have the advantage such as being flexible and easy to apply on nonuniform grids. We introduce an ANN based troubled-cell indicator by constructing a multilayer perceptron (MLP) model, one of the most common ANN models. The third-order WENO scheme is focused in this paper. Extensive numerical experiments for solving various scalar equations with both convex and non-convex cases, and the Euler systems of equations on uniform and nonuniform grids of one-dimensional (1D) and two-dimensional (2D) domains, are performed to show the accuracy and nonlinear stability of the proposed hybrid finite volume WENO scheme with the MLP troubled-cell indicator. Significant accuracy improvement and computational-cost saving over the original WENO scheme are observed. Numerical experiments and comparisons with the widely-used KXRCF indicator also show the good performance of the MLP troubled-cell indicator. Although the MLP troubled-cell indicator is trained on uniform grids, it performs very well on nonuniform grids obtained by randomly perturbing uniform grids.

Abstract. The Lagrangian method has attracted increasing interest in multimaterial fluid flow fields, as it allows the grid to move with the fluid velocity, and it can track the material interfaces automatically and sharply. There are many studies discussing the finite volume (FV) or discontinuous Galerkin (DG) Lagrangian schemes, but little attention has been paid to the finite difference Lagrangian schemes. In this paper, we construct a high-order conservative finite difference alternative weighted essentially nonoscillatory (AWENO) Lagrangian scheme for multidimensional compressible Euler equations in curvilinear coordinates. This is the first high-order finite difference work in pure Lagrangian framework. Specifically, the governing equations are represented in Lagrangian coordinates and are truly hyperbolic. Then we design a fifth-order finite difference AWENO scheme with third-order strong stability preserving Runge–Kutta time-marching. This scheme is conservative, and it can not only achieve arbitrary high-order accuracy in smooth regions, but also avoid numerical oscillations near discontinuities. Moreover, by appropriately discretizing the metric derivatives, the discrete compatibility conditions and free-stream preserving property are proved. More importantly, there is no mass exchange across the moving grid points. Compared to the FV or DG Lagrangian schemes, our finite difference Lagrangian scheme is simpler in coding, it easier to achieve arbitrary high-order accuracy, and, most importantly, it is easily extended and is more efficient in high dimensions. A series of numerical experiments including one-dimensional, two-dimensional, and three-dimensional examples are given to demonstrate the good performance of the scheme in terms of high-order accuracy, nonoscillation, high resolution for discontinuities, and free-stream preservation.

We combine Patankar-type methods with suitable relaxation procedures that are capable of ensuring correct dissipation or conservation of functionals such as entropy or energy while producing unconditionally positive and conservative approximations. To that end, we adapt the relaxation algorithm to enforce positivity by using either ideas from the dense output framework when a linear invariant must be preserved, or simply a geometric mean if the only constraint is positivity preservation. The latter merely requires the solution of a scalar nonlinear equation while former results in a coupled linear-nonlinear system of equations. We present sufficient conditions for the solvability of the respective equations. Several applications in the context of ordinary and partial differential equations are presented, and the theoretical findings are validated numerically.

We prove the stability and convergence of the high order discontinuous Galerkin scheme to spherically symmetric Einstein-scalar equations for a class of large initial data that ensures the formation of a black hole. Having chosen the Bondi coordinate system, we achieve L2 stability and obtain the optimal error estimates.

We present a struture-preserving solver for particle-wave interaction in magnetized plasmas. The solver combines a conservative local discontinuous Galerkin (LDG) scheme for the interaction part with a trajectory averaging method for the Hamiltonian flow part. The proposed LDG scheme is an extension of the conservative scheme we developed in 2023. The trajectory averaging method significantly reduces computational cost by taking advantage of the multiscale feature of this system. By introducing a novel concept “trajectory bundle”, we transform a continuous topological problem into a discrete graph theory problem. Numerical examples for a non-uniform magnetized plasma in an infinitely long symmetric cylinder is presented. It is verified that the connection-proportion algorithm allows to distinguish different trajectory bundles, and the proposed DG scheme rigorously preserves all the conservation laws.

Numerically solving magnetohydrodynamic (MHD) equations faces many challenges: avoiding divergence error, maintaining positivity, and satisfying entropy conditions. Among discontinuous Galerkin (DG) schemes, there has been a modal version that is locally divergence-free and positivity preserving and a nodal version that is semi-discretely entropy stable. In this work, we develop a DG scheme that combines the advantages of these two and solves all the three challenges. The key ingredients that bring these two schemes together are an HLL numerical flux with entropy stable signal speed estimates and a locally divergence-free projection. To handle problems with strong shocks, the essentially oscillation-free damping is applied. Various numerical experiments verify the accuracy and robustness of our method.

A local discontinuous Galerkin (LDG) scheme coupled with a third-order strong stability-preserving (SSP) Runge–Kutta time integration, is proposed for the numerical simulation of non-Fourier heat transfer in longitudinal fins subject to temperature-dependent convective and radiative heat transfer incorporating direct interaction between the fin surface and the base. The mathematical model is based on the Maxwell–Cattaneo–Vernotte (MCV) hyperbolic heat conduction framework, which accounts for the finite speed of thermal wave propagation via a relaxation time parameter. Three representative fin profiles — trapezoidal, rectangular, and dovetail — are considered by varying the taper ratio. The scheme is verified against exact solutions and achieves third-order accuracy. Moreover, it accurately captures the sharp wavefront without oscillations. The effects of the key physical parameters on the transient temperature distribution and fin efficiency are systematically examined. It is found that the thermal wave speed is primarily governed by the relaxation parameter Ve and the thermal conductivity parameter β, but is independent of the surface heat transfer mechanisms and fin geometry. Both convection and radiation effects play an important role in the process of heat transfer, and strengthening these effects will enhance heat loss and increase heat transfer rate. The time-averaged fin efficiency increases with Ve and is consistently higher for the dovetail fin than for both the rectangular and trapezoidal fins under the same operating conditions.

Solutions of hyperbolic conservation laws exhibit both smooth structures across large scales and sharp localized features such as shocks and contact discontinuities, making them difficult to approximate accurately with existing neural operators. The Fourier Neural Operator (FNO) captures long-range interactions well but tends to smear localized structures through excessive numerical dissipation. To address this, we propose a Local-Global Neural Operator (LGNO) that learns a one-step discrete flow map by combining a global FNO branch for representing smooth dynamics at large scales with a local multiresolution branch for enhancing localized discontinuities and nonsmooth features. The model is trained with a one-step loss that combines a physical space prediction term and a spectral penalty on high frequencies to suppress spurious oscillations near steep fronts. On a large collection of benchmarks in one and two dimensions, LGNO consistently outperforms FNO baselines with matched parameter counts, reducing one-step errors by factors of 2-5 and remaining significantly more accurate over long autoregressive rollouts. Most strikingly, although it is trained only on short-time data from a high-order WENO-Z scheme, the long-time rollout of LGNO on a coarse $256^2$ grid exhibits lower numerical dissipation than the same WENO-Z scheme run on a finer $512^2$ grid, while being orders of magnitude cheaper to evaluate. These results suggest that, with an appropriate architecture and training objective, learned operators can effectively learn discrete flow maps. They further suggest that such learned operators have the potential to control long-time numerical dissipation better than the conventional shock-capturing schemes that generate the training data.

Entropy inequalities are fundamental to the well-posedness of hyperbolic conservation laws, providing the essential criterion for selecting the physically admissible solution among infinitely many weak solutions. Chen and Shu [J. Comput. Phys. 345 (2017)] proposed a unified framework for constructing high-order discontinuous Galerkin (DG) methods that satisfy entropy inequalities for any given entropy via specific numerical quadrature; however, their accompanying error analysis was limited to the truncation error level, leaving a critical gap in the rigorous convergence theory for these entropy-stable schemes. This paper closes that gap by establishing rigorous a priori error estimates for semi-discrete entropy-stable DG (ESDG) methods on general unstructured meshes for hyperbolic conservation laws. The analysis applies to both scalar equations and systems, and is built upon a finite-difference-type consistency-stability argument carried out directly at the nodal level. Under a polynomial-reconstruction hypothesis and an L^∞ a priori bound, we prove an O(h^k) error estimate in a quadrature-based norm, which is equivalent to the broken L^2 norm on the finite-dimensional reconstruction space. We further extend this framework to the entropy-stable oscillation-free DG (ESOFDG) method introduced by Liu, Lu, and Shu [SIAM J. Sci. Comput. 46 (2024)], demonstrating that the additional damping terms do not degrade the convergence order. Numerical experiments suggest that the observed convergence rates may exceed the theoretical prediction by up to half an order.

We construct a non-polynomial local discontinuous Galerkin (LDG) scheme for the prescribed mean curvature equation to approximate boundary gradient blow-up solutions and obtain error estimates.

In this paper, a stable and arbitrary high order spectral volume (SV) method is proposed for solving linear hyperbolic conservation laws on unstructured quadrilateral meshes. The SV scheme is constructed with subdivision points being zeros of a class of parameterized polynomials. A unified proof for the $L^2$ stability of the SV method is provided under conditions that the parameter $c>-\frac{1}{k(k+1)}$ and the underlying mesh is an $h^{1+\gamma}, \gamma\ge 1$ parallelogram mesh. The optimal a priori error estimate of SV methods is established. Numerical experiments are presented to verify all theoretical findings.

We establish a simple, rigorous, and easy to implement connection between the classical continuous finite element method (FEM) and the discontinuous Galerkin (DG) method for Poisson's problem. The key idea is to insert a vanishing-thickness layer of "dummy" elements along cell interfaces. By modifying the diffusion coefficient on these elements to be proportional to their thickness, we prove the FEM formulation converges to Babuška-Zlámal DG with trapezoidal edge quadrature. The scheme is trivial to implement by (i) a mesh edit that introduces degenerate interface elements and (ii) a single Jacobian threshold in an otherwise unmodified FEM code to handle the degenerate elements via the tempered finite element (TFEM) framework. We provide a rigorous derivation of the resulting TFEM-DG scheme, prove optimal H^1 and L^2 error estimates, and present numerical experiments in 2D and 3D. The method allows for simple implementation of DG in a FEM code and even adaptive element-by-element switching between FEM and DG with minimal coding effort. The framework is readily extensible, as we will demonstrate in a companion paper dedicated to evolutionary nonlinear first-order hyperbolic systems.

Stochastic phenomena are commonly observed in pedestrian flow. However, the existing models for pedestrian dynamics rely on averaged inputs and yield deterministic outputs only, and thereby fail to capture the stochastic variabilities inherent in pedestrian dynamics. This study builds upon Hughes’ dynamic continuum model to develop mathematical models for stochastic pedestrian dynamics that explicitly consider two types of stochastic characteristics: demand stochasticity and behavioral stochasticity. The proposed system, represented as a set of time-dependent stochastic partial differential equations, is solved using a combination of the Monte Carlo (MC) method or Quasi Monte Carlo (QMC) method and efficient numerical schemes, such as a fifth-order weighted essentially non-oscillatory finite difference scheme and the fast sweeping method. Benchmarking scenarios are designed and simulated, and the numerical results demonstrate the convergence and computational performance of the MC and QMC methods. The advantage of stochastic modeling is evident given the significant differences between the averaged stochastic outputs and deterministic outputs, attributable to the strength of stochasticity and extent of stochastic dimensions. Moreover, based on stochastic data inputs, the proposed stochastic models and numerical solutions can clarify the probabilistic distributions of key indicators, such as density, which are valuable for the design and improvement of pedestrian facilities.

To address the computational bottlenecks in efficient, accurate, and robust simulations of the Richards' equation for unsaturated flow in geotechnical media, we propose a third-order explicit-implicit-null (EIN) time integration method. For spatial discretization, a conservative high-order finite-difference multi-resolution weighted essentially non-oscillatory (WENO) scheme is constructed to accurately resolve sharp wetting fronts and suppress spurious oscillations. The core idea of the EIN framework is to add and subtract a suitably large Laplacian operator at one side of the Richards' equation, followed by an implicit-explicit (IMEX) time-marching strategy that treats the added linear diffusion term implicitly and all other terms explicitly. This treatment significantly relaxes the time-step restriction, avoids costly nonlinear iterations, and thus substantially improves computational efficiency. To further ensure physical rationality and numerical robustness, a bound-preserving sweeping technique is incorporated to maintain solutions within physical bounds, while preserving conservation and accuracy. The proposed algorithm is easy to implement and readily integrable into geotechnical applications. Numerical experiments on canonical infiltration benchmarks and challenging three-dimensional (3D) application problems-including layered and heterogeneous soils, and a realistic case with field-sampled black soil parameters from Northeast China-demonstrate the performance and practical applicability, achieving over 90% computational savings compared to the explicit Runge-Kutta method for large-scale and 3D simulations.

The diffusive-viscous wave equations (DVWE) arise in geophysics and describe the propagation of seismic waves in fluid-saturated media. In this paper, we investigate three implicit-explicit Runge-Kutta time discretization schemes up to third-order, coupled with local discontinuous Galerkin methods, for solving the DVWE with variable coefficients. The unconditional energy stability of the proposed schemes is rigorously established through energy analysis. The main technique exploits the relationship between auxiliary and primary variables in the general case where the auxiliary variable is defined as the gradient of the primary variable multiplied by a variable coefficient. Moreover, by utilizing elliptic projection together with the supercloseness property, we successfully achieve optimal error estimates in the L2 norm for the DVWE with general smooth variable coefficients. Numerical experiments validate the theoretical results and demonstrate the performance of the proposed schemes for both smooth and discontinuous media.

This paper develops a local discontinuous Galerkin (LDG) method with semi-implicit time discretization for non-Fourier heat transfer governed by the Maxwell-Cattaneo-Vernotte (MCV) model in longitudinal fins with arbitrary profiles (1-x)p (p = 0 to 3), subjected to both cosinusoidal and periodic square wave thermal excitations. The method demonstrates excellent discontinuity-capturing capability and robust performance for stiff systems with large values of p. The study reveals that wave speed in the MCV framework is determined solely by the dimensionless relaxation parameter beta, exhibiting no dependence on fin profile or excitation characteristics. Among the profiles examined, rectangular fins (p = 0) exhibit the highest thermal efficiency, which decreases monotonically with increasing p. Between the two excitation modes, square excitation yields slightly superior performance compared to cosinusoidal wave excitation. Moreover, larger relaxation parameters may enhance time-averaged efficiency under sustained periodic heating but at the cost of slower transient response. These findings provide fundamental insights for designing thermal management systems under hyperbolic heat conduction.

Accurate and efficient reconstruction techniques are essential in multiresolution analysis and image compression, particularly when the data are represented as cell averages. In this work, we present a non-separable progressive multivariate Weighted Essentially Non-Oscillatory (WENO) scheme specifically designed for cell-average data, with applications to digital image processing. The proposed method extends Harten’s multiresolution framework through a non-linear WENO reconstruction adapted to the cell-average context. In contrast to classical WENO schemes, the progressive strategy allows the recursive recovery of high-order accuracy even when the largest stencil is affected by a discontinuity. The method achieves high-order accuracy in smooth regions together with stable, non-oscillatory behavior near discontinuities. We also establish theoretical results regarding the consistency and approximation properties of the method. Finally, several numerical experiments on piecewise smooth functions and digital images are presented to demonstrate its performance and validate its effectiveness against the linear Lagrange reconstruction of the same order of accuracy.

This paper presents a novel two-dimensional intersection-based remapping method for isoparametric curvilinear meshes within the indirect arbitrary Lagrangian-Eulerian (ALE) framework, addressing the challenges of transferring physical quantities between high-order curved-edge meshes. Our method leverages the Weiler-Atherton clipping algorithm to compute intersections between curved-edge quadrangles, enabling robust handling of arbitrary order isoparametric curves. By integrating multi-resolution weighted essentially non-oscillatory (WENO) reconstruction, we achieve high-order accuracy while suppressing numerical oscillations near discontinuities. A positivity-preserving limiter is further applied to ensure physical quantities such as density remain non-negative without compromising conservation or accuracy. Notably, the computational cost of handling higher-order curved meshes, such as cubic or even higher-degree parametric curves, does not significantly increase compared to second-order curved meshes. This ensures that our method remains efficient and scalable, making it applicable to arbitrary two-dimensional high-order isoparametric curvilinear cells without compromising performance. Numerical experiments demonstrate that the proposed method achieves high-order accuracy, strict conservation (with errors approaching machine precision), essential non-oscillation and positivity-preserving. The proposed approach is currently restricted to two-dimensional meshes, and an extension to fully three-dimensional curved polyhedral mesh is beyond the scope of the present work.