Kailiang Wu

Kailiang Wu, Chi-Wang Shu

This paper proposes and analyzes arbitrarily high-order discontinuous Galerkin (DG) and finite volume methods which provably preserve the positivity of density and pressure for the ideal MHD on general meshes. Unified auxiliary theories are built for rigorously analyzing the positivity-preserving (PP) property of MHD schemes with a HLL type flux on polytopal meshes in any space dimension. The main challenges overcome here include establishing relation between the PP property and discrete divergence of magnetic field on general meshes, and estimating proper wave speeds in the HLL flux to ensure the PP property. In 1D case, we prove that the standard DG and finite volume methods with the proposed HLL flux are PP, under condition accessible by a PP limiter. For multidimensional conservative MHD system, standard DG methods with a PP limiter are not PP in general, due to the effect of unavoidable divergence-error. We construct provably PP high-order DG and finite volume schemes by proper discretization of symmetrizable MHD system, with two divergence-controlling techniques: locally divergence-free elements and a penalty term. The former leads to zero divergence within each cell, while the latter controls the divergence error across cell interfaces. Our analysis reveals that a coupling of them is important for positivity preservation, as they exactly contribute the discrete divergence-terms absent in standard DG schemes but crucial for ensuring the PP property. Numerical tests confirm the PP property and the effectiveness of proposed PP schemes. Unlike conservative MHD system, the exact smooth solutions of symmetrizable MHD system are proved to retain the positivity even if the divergence-free condition is not satisfied. Our analysis and findings further the understanding, at both discrete and continuous levels, of the relation between the PP property and the divergence-free constraint.

Kailiang Wu, Chi-Wang Shu

The density and pressure are positive physical quantities in magnetohydrodynamics (MHD). Design of provably positivity-preserving (PP) numerical schemes for ideal compressible MHD is highly desirable, but remains a challenge especially in the multidimensional cases. In this paper, we first develop uniformly high-order discontinuous Galerkin (DG) schemes which provably preserve the positivity of density and pressure for multidimensional ideal MHD. The schemes are constructed by using the locally divergence-free DG schemes for the symmetrizable ideal MHD equations as the base schemes, a PP limiter to enforce the positivity of the DG solutions, and the strong stability preserving methods for time discretization. The significant innovation is that we discover and rigorously prove the PP property of the proposed DG schemes by using a novel equivalent form of the admissible state set and very technical estimates. Several two-dimensional numerical examples further confirm the PP property, and demonstrate the accuracy, effectiveness and robustness of the proposed PP methods.

Ce Zhang, Kailiang Wu, Zhihai He · cmu

Given an unknown dynamical system, what is the minimum number of samples needed for effective learning of its governing laws and accurate prediction of its future evolution behavior, and how to select these critical samples? In this work, we propose to explore this problem based on a design approach. Starting from a small initial set of samples, we adaptively discover critical samples to achieve increasingly accurate learning of the system evolution. One central challenge here is that we do not know the network modeling error since the ground-truth system state is unknown, which is however needed for critical sampling. To address this challenge, we introduce a multi-step reciprocal prediction network where forward and backward evolution networks are designed to learn the temporal evolution behavior in the forward and backward time directions, respectively. Very interestingly, we find that the desired network modeling error is highly correlated with the multi-step reciprocal prediction error, which can be directly computed from the current system state. This allows us to perform a dynamic selection of critical samples from regions with high network modeling errors for dynamical systems. Additionally, a joint spatial-temporal evolution network is introduced which incorporates spatial dynamics modeling into the temporal evolution prediction for robust learning of the system evolution operator with few samples. Our extensive experimental results demonstrate that our proposed method is able to dramatically reduce the number of samples needed for effective learning and accurate prediction of evolution behaviors of unknown dynamical systems by up to hundreds of times.

Kailiang Wu, Huazhong Tang, Dongbin Xiu

This paper is concerned with generalized polynomial chaos (gPC) approximation for a general system of quasilinear hyperbolic conservation laws with uncertainty. The one-dimensional (1D) hyperbolic system is first symmetrized with the aid of left eigenvector matrix of the Jacobian matrix. Stochastic Galerkin method is then applied to derive the equations for the gPC expansion coefficients. The resulting deterministic gPC Galerkin system is proved to be symmetrically hyperbolic. This important property then allows one to use a variety of numerical schemes for spatial and temporal discretization. Here a higher-order and path-conservative finite volume WENO scheme is adopted in space, along with a third-order total variation diminishing Runge-Kutta method in time. The method is further extended to two-dimensional (2D) quasilinear hyperbolic system with uncertainty, where the symmetric hyperbolicity of the one-dimensional system is carried over via the operator splitting technique. Several 1D and 2D numerical experiments are conducted to demonstrate the accuracy and effectiveness of the proposed gPC stochastic Galerkin method.

Kailiang Wu, Huazhong Tang

This paper first studies the admissible state set $\mathcal G$ of relativistic magnetohydrodynamics (RMHD). It paves a way for developing physical-constraints-preserving (PCP) schemes for RMHD equations with the solutions in $\mathcal G$. To overcome the difficulties arising from the extremely strong nonlinearities and no explicit formulas of the primitive variables and flux vectors with respect to the conservative vector, two equivalent forms of $\mathcal G$ with explicit constraints on the conservative vector are skillfully discovered. The first is derived by analyzing roots of several polynomials and transferring successively them, and further used to prove the convexity of $\mathcal G$ with the aid of semi-positive definiteness of the second fundamental form of a hypersurface. While the second is derived based on the convexity and then used to show the orthogonal invariance of $\mathcal G$. The Lax-Friedrichs (LxF) splitting property does not hold generally for nonzero magnetic field, but by a constructive inequality and pivotal techniques, we discover the generalized LxF splitting properties, combining the convex combination of some LxF splitting terms with a discrete divergence-free condition of magnetic field. Based on above analyses, several one- and two-dimensional PCP schemes are then studied. In 1D case, a first-order accurate LxF type scheme is first proved to be PCP under the CFL condition, and then high-order PCP schemes are proposed via a PCP limiter. In 2D case, the discrete divergence-free condition and PCP property are analyzed for a first-order LxF scheme, and two sufficient conditions are derived for high-order accurate PCP schemes. Our analysis reveals in theory for the first time that the discrete divergence-free condition is closely connected with the PCP property. Several numerical examples demonstrate theoretical findings and performance of numerical schemes.

Kailiang Wu

Numerical schemes provably preserving the positivity of density and pressure are highly desirable for MHD, but the rigorous positivity-preserving (PP) analysis remains challenging. The difficulties mainly arise from the intrinsic complexity of the MHD equations as well as the indeterminate relation between the PP property and the divergence-free condition on magnetic field. We present the first rigorous PP analysis of conservative schemes with Lax-Friedrichs (LF) flux for ideal MHD. The significant innovation is the discovery of theoretical connection between PP property and a discrete divergence-free (DDF) condition. This connection is established through the generalized LF splitting properties, which are alternatives of the usually-expected LF splitting property that does not hold for ideal MHD. The generalized LF splitting properties involve a number of admissible states strongly coupled by DDF condition, making their derivation very difficult. We derive these properties via a novel equivalent form of the admissible state set and an important inequality skillfully constructed by technical estimates. Rigorous PP analysis is presented for finite volume and discontinuous Galerkin schemes with LF flux on uniform Cartesian meshes. In 1D case, PP property is proved for the first-order scheme with proper numerical viscosity, and also for arbitrarily high-order schemes under conditions accessible by a PP limiter. In 2D case, we show that the DDF condition is necessary and crucial for achieving PP property. It is observed that even slightly violating the proposed DDF condition may cause failure to preserve the positivity of pressure. We prove that the 2D LF type scheme with proper numerical viscosity preserves the positivity and DDF condition. Sufficient conditions are derived for 2D PP high-order schemes, and extension to 3D is discussed. Numerical examples confirm the theoretical findings.

Kailiang Wu, Huazhong Tang

The paper proposes a second-order accurate direct Eulerian generalized Riemann problem (GRP) scheme for the spherically symmetric general relativistic hydrodynamical (RHD) equations and a second-order accurate discretization for the spherically symmetric Einstein (SSE) equations. The former is directly using the Riemann invariants and the Runkine-Hugoniot jump conditions to analytically resolve the left and right nonlinear waves of the local GRP in the Eulerian formulation together with the local change of the metrics to obtain the limiting values of the time derivatives of the conservative variables along the cell interface and the numerical flux for the GRP scheme. While the latter utilizes the energy-momentum tensor obtained in the GRP solver to evaluate the fluid variables in the SSE equations and keeps the continuity of the metrics at the cell interfaces. Several numerical experiments show that the GRP scheme can achieve second-order accuracy and high resolution, and is effective for spherically symmetric general RHD problems.

Kailiang Wu

The paper develops high-order physical-constraint-preserving (PCP) methods for general relativistic hydrodynamic (GRHD) equations, equipped with a general equation of state. Here the physical constraints, describing the admissible states of GRHD, are referred to the subluminal constraint on the fluid velocity and the positivity of the density, pressure and specific internal energy. Preserving these constraints is very important for robust computations, otherwise violating one of them will lead to the ill-posed problem and numerical instability. To overcome the difficulties arising from the inherent strong nonlinearity contained in the constraints, we derive an equivalent definition of the admissible states. Using this definition, we prove the convexity, scaling invariance and Lax-Friedrichs (LxF) splitting property of the admissible state set $\mathcal G$, and discover the dependence of $\mathcal G$ on the spacetime metric. Unfortunately, such dependence yields the non-equivalence of $\mathcal G$ at different points in curved spacetime, and invalidates the convexity of $\mathcal G$ in analyzing PCP schemes. This obstacle is effectively overcame by introducing a new formulation of the GRHD equations. Based on this formulation and the above theories, a first-order LxF scheme is designed on general unstructured mesh and rigorously proved to be PCP under a CFL condition. With two types of PCP limiting procedures, we design high-order, {\em provably} (not probably) PCP methods under discretization on the proposed new formulation. These high-order methods include the PCP finite difference, finite volume and discontinuous Galerkin methods.

Kailiang Wu, Huazhong Tang

The paper studies the physical-constraints-preserving (PCP) schemes for multi-dimensional special relativistic magnetohydrodynamics with a general equation of state (EOS) on more general meshes. It is an extension of the work [Math. Models Methods Appl. Sci., 27:1871-1928, 2017] which focuses on the ideal EOS and uniform Cartesian meshes. The general EOS without a special expression poses some additional difficulties in discussing the mathematical properties of admissible state set with the physical constraints on the fluid velocity, density and pressure. Rigorous analyses are provided for the PCP property of finite volume or discontinuous Galerkin schemes with the Lax-Friedrichs (LxF) type flux on a general mesh with non-self-intersecting polytopes. Those are built on a more general form of generalized LxF splitting property and a different convex decomposition technique. It is shown in theory that the PCP property is closely connected with a discrete divergence-free condition, which is proposed on the general mesh and milder than that in [Math. Models Methods Appl. Sci., 27:1871-1928, 2017].

Junfeng Chen, Kailiang Wu

This paper proposes a novel deep learning approach for learning operators in semigroup, with applications to modeling unknown autonomous dynamical systems using time series data collected at varied time lags. It is a sequel to the previous flow map learning (FML) works [T. Qin, K. Wu, and D. Xiu, J. Comput. Phys., 395:620--635, 2019], [K. Wu and D. Xiu, J. Comput. Phys., 408:109307, 2020], and [Z. Chen, V. Churchill, K. Wu, and D. Xiu, J. Comput. Phys., 449:110782, 2022], which focused on learning single evolution operator with a fixed time step. This paper aims to learn a family of evolution operators with variable time steps, which constitute a semigroup for an autonomous system. The semigroup property is very crucial and links the system's evolutionary behaviors across varying time scales, but it was not considered in the previous works. We propose for the first time a framework of embedding the semigroup property into the data-driven learning process, through a novel neural network architecture and new loss functions. The framework is very feasible, can be combined with any suitable neural networks, and is applicable to learning general autonomous ODEs and PDEs. We present the rigorous error estimates and variance analysis to understand the prediction accuracy and robustness of our approach, showing the remarkable advantages of semigroup awareness in our model. Moreover, our approach allows one to arbitrarily choose the time steps for prediction and ensures that the predicted results are well self-matched and consistent. Extensive numerical experiments demonstrate that embedding the semigroup property notably reduces the data dependency of deep learning models and greatly improves the accuracy, robustness, and stability for long-time prediction.

Junfeng Chen, Kailiang Wu, Dongbin Xiu

Equations, particularly differential equations, are fundamental for understanding natural phenomena and predicting complex dynamics across various scientific and engineering disciplines. However, the governing equations for many complex systems remain unknown due to intricate underlying mechanisms. Recent advancements in machine learning and data science offer a new paradigm for modeling unknown equations from measurement or simulation data. This paradigm shift, known as data-driven discovery or modeling, stands at the forefront of AI for science, with significant progress made in recent years. In this paper, we introduce a systematic framework for data-driven modeling of unknown equations using deep learning. This versatile framework is capable of learning unknown ODEs, PDEs, DAEs, IDEs, SDEs, reduced or partially observed systems, and non-autonomous differential equations. Based on this framework, we have developed Deep Unknown Equations (DUE), an open-source software package designed to facilitate the data-driven modeling of unknown equations using modern deep learning techniques. DUE serves as an educational tool for classroom instruction, enabling students and newcomers to gain hands-on experience with differential equations, data-driven modeling, and contemporary deep learning approaches such as FNN, ResNet, generalized ResNet, operator semigroup networks (OSG-Net), and Transformers. Additionally, DUE is a versatile and accessible toolkit for researchers across various scientific and engineering fields. It is applicable not only for learning unknown equations from data but also for surrogate modeling of known, yet complex, equations that are costly to solve using traditional numerical methods. We provide detailed descriptions of DUE and demonstrate its capabilities through diverse examples, which serve as templates that can be easily adapted for other applications.

Junfeng Chen, Kailiang Wu

Operator learning for Partial Differential Equations (PDEs) is rapidly emerging as a promising approach for surrogate modeling of intricate systems. Transformers with the self-attention mechanism$\unicode{x2013}$a powerful tool originally designed for natural language processing$\unicode{x2013}$have recently been adapted for operator learning. However, they confront challenges, including high computational demands and limited interpretability. This raises a critical question: Is there a more efficient attention mechanism for Transformer-based operator learning? This paper proposes the Position-induced Transformer (PiT), built on an innovative position-attention mechanism, which demonstrates significant advantages over the classical self-attention in operator learning. Position-attention draws inspiration from numerical methods for PDEs. Different from self-attention, position-attention is induced by only the spatial interrelations of sampling positions for input functions of the operators, and does not rely on the input function values themselves, thereby greatly boosting efficiency. PiT exhibits superior performance over current state-of-the-art neural operators in a variety of complex operator learning tasks across diverse PDE benchmarks. Additionally, PiT possesses an enhanced discretization convergence feature, compared to the widely-used Fourier neural operator.

Kailiang Wu

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 {\em 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. {\em 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.

Zhen Chen, Victor Churchill, Kailiang Wu et al.

We present a numerical framework for deep neural network (DNN) modeling of unknown time-dependent partial differential equations (PDE) using their trajectory data. Unlike the recent work of [Wu and Xiu, J. Comput. Phys. 2020], where the learning takes place in modal/Fourier space, the current method conducts the learning and modeling in physical space and uses measurement data as nodal values. We present a DNN structure that has a direct correspondence to the evolution operator of the underlying PDE, thus establishing the existence of the DNN model. The DNN model also does not require any geometric information of the data nodes. Consequently, a trained DNN defines a predictive model for the underlying unknown PDE over structureless grids. A set of examples, including linear and nonlinear scalar PDE, system of PDEs, in both one dimension and two dimensions, over structured and unstructured grids, are presented to demonstrate the effectiveness of the proposed DNN modeling. Extension to other equations such as differential-integral equations is also discussed.

Zhen Chen, Kailiang Wu, Dongbin Xiu

We study the problem of identifying unknown processes embedded in time-dependent partial differential equation (PDE) using observational data, with an application to advection-diffusion type PDE. We first conduct theoretical analysis and derive conditions to ensure the solvability of the problem. We then present a set of numerical approaches, including Galerkin type algorithm and collocation type algorithm. Analysis of the algorithms are presented, along with their implementation detail. The Galerkin algorithm is more suitable for practical situations, particularly those with noisy data, as it avoids using derivative/gradient data. Various numerical examples are then presented to demonstrate the performance and properties of the numerical methods.

Jun Hou, Tong Qin, Kailiang Wu et al.

A novel correction algorithm is proposed for multi-class classification problems with corrupted training data. The algorithm is non-intrusive, in the sense that it post-processes a trained classification model by adding a correction procedure to the model prediction. The correction procedure can be coupled with any approximators, such as logistic regression, neural networks of various architectures, etc. When training dataset is sufficiently large, we prove that the corrected models deliver correct classification results as if there is no corruption in the training data. For datasets of finite size, the corrected models produce significantly better recovery results, compared to the models without the correction algorithm. All of the theoretical findings in the paper are verified by our numerical examples.

Kailiang Wu, Dongbin Xiu

We present a framework for recovering/approximating unknown time-dependent partial differential equation (PDE) using its solution data. Instead of identifying the terms in the underlying PDE, we seek to approximate the evolution operator of the underlying PDE numerically. The evolution operator of the PDE, defined in infinite-dimensional space, maps the solution from a current time to a future time and completely characterizes the solution evolution of the underlying unknown PDE. Our recovery strategy relies on approximation of the evolution operator in a properly defined modal space, i.e., generalized Fourier space, in order to reduce the problem to finite dimensions. The finite dimensional approximation is then accomplished by training a deep neural network structure, which is based on residual network (ResNet), using the given data. Error analysis is provided to illustrate the predictive accuracy of the proposed method. A set of examples of different types of PDEs, including inviscid Burgers' equation that develops discontinuity in its solution, are presented to demonstrate the effectiveness of the proposed method.

Kailiang Wu, Tong Qin, Dongbin Xiu

We present a numerical approach for approximating unknown Hamiltonian systems using observation data. A distinct feature of the proposed method is that it is structure-preserving, in the sense that it enforces conservation of the reconstructed Hamiltonian. This is achieved by directly approximating the underlying unknown Hamiltonian, rather than the right-hand-side of the governing equations. We present the technical details of the proposed algorithm and its error estimate in a special case, along with a practical de-noising procedure to cope with noisy data. A set of numerical examples are then presented to demonstrate the structure-preserving property and effectiveness of the algorithm.

Tong Qin, Kailiang Wu, Dongbin Xiu

We present a numerical framework for approximating unknown governing equations using observation data and deep neural networks (DNN). In particular, we propose to use residual network (ResNet) as the basic building block for equation approximation. We demonstrate that the ResNet block can be considered as a one-step method that is exact in temporal integration. We then present two multi-step methods, recurrent ResNet (RT-ResNet) method and recursive ReNet (RS-ResNet) method. The RT-ResNet is a multi-step method on uniform time steps, whereas the RS-ResNet is an adaptive multi-step method using variable time steps. All three methods presented here are based on integral form of the underlying dynamical system. As a result, they do not require time derivative data for equation recovery and can cope with relatively coarsely distributed trajectory data. Several numerical examples are presented to demonstrate the performance of the methods.

Kailiang Wu, Dongbin Xiu

We present effective numerical algorithms for locally recovering unknown governing differential equations from measurement data. We employ a set of standard basis functions, e.g., polynomials, to approximate the governing equation with high accuracy. Upon recasting the problem into a function approximation problem, we discuss several important aspects for accurate approximation. Most notably, we discuss the importance of using a large number of short bursts of trajectory data, rather than using data from a single long trajectory. Several options for the numerical algorithms to perform accurate approximation are then presented, along with an error estimate of the final equation approximation. We then present an extensive set of numerical examples of both linear and nonlinear systems to demonstrate the properties and effectiveness of our equation recovery algorithms.

Kailiang Wu, Dongbin Xiu

We present an explicit construction for feedforward neural network (FNN), which provides a piecewise constant approximation for multivariate functions. The proposed FNN has two hidden layers, where the weights and thresholds are explicitly defined and do not require numerical optimization for training. Unlike most of the existing work on explicit FNN construction, the proposed FNN does not rely on tensor structure in multiple dimensions. Instead, it automatically creates Voronoi tessellation of the domain, based on the given data of the target function, and piecewise constant approximation of the function. This makes the construction more practical for applications. We present both theoretical analysis and numerical examples to demonstrate its properties.

Yangyu Kuang, Kailiang Wu, Huazhong Tang

The paper develops high order accurate Runge-Kutta discontinuous local evolution Galerkin (RKDLEG) methods on the cubed-sphere grid for the shallow water equations (SWEs). Instead of using the dimensional splitting method or solving one-dimensional Riemann problem in the direction normal to the cell interface, the RKDLEG methods are built on genuinely multi-dimensional approximate local evolution operator of the locally linearized SWEs on a sphere by considering all bicharacteristic directions. Several numerical experiments are conducted to demonstrate the accuracy and performance of our RKDLEG methods, in comparison to the Runge-Kutta discontinuous Galerkin method with Godunov's flux etc.

Kailiang Wu, Huazhong Tang

The paper develops a Newton multigrid (MG) method for one- and two-dimensional steady-state shallow water equations (SWEs) with topography and dry areas.It solves the nonlinear system arising from the well-balanced finite volume discretization of the steady-state SWEs by using Newton's method as the outer iteration and a geometric MG method with the block symmetric Gauss-Seidel smoother as the inner iteration. The proposed Newton MG method makes use of the local residual to regularize the Jacobian matrix of the Newton iteration, and can handle the steady-state problem with wet/dry transitions. Several numerical experiments are conducted to demonstrate the efficiency, robustness, and well-balanced property of the proposed method. The relation between the convergence behavior of the Newton MG method and the distribution of the eigenvalues of the iteration matrix is detailedly discussed.

Kailiang Wu, Huazhong Tang

The ideal gas equation of state (EOS) with a constant adiabatic index is a poor approximation for most relativistic astrophysical flows, although it is commonly used in relativistic hydrodynamics. The paper develops high-order accurate physical-constraints-preserving (PCP) central discontinuous Galerkin (DG) methods for the one- and two-dimensional special relativistic hydrodynamic (RHD) equations with a general EOS. It is built on the theoretical analysis of the admissible states for the RHD and the PCP limiting procedure enforcing the admissibility of central DG solutions. The convexity, scaling and orthogonal invariance, and Lax-Friedrichs splitting property of the admissible state set are first proved with the aid of its equivalent form, and then the high-order central DG methods with the PCP limiting procedure and strong stability preserving time discretization are proved to preserve the positivity of the density, pressure, and specific internal energy, and the bound of the fluid velocity, maintain the high-order accuracy, and be $L^1$-stable. The accuracy, robustness, and effectiveness of the proposed methods are demonstrated by several 1D and 2D numerical examples involving large Lorentz factor, strong discontinuities, or low density or pressure etc.