Search

In this paper, we study the uniform accuracy of implicit-explicit (IMEX) Runge-Kutta (RK) schemes for general linear hyperbolic relaxation systems satisfying the structural stability condition proposed in \cite{yong_singular_1999}. We establish the uniform stability and accuracy of a class of IMEX-RK schemes with spatial discretization using a Fourier spectral method. Our results demonstrate that the accuracy of the fully discretized schemes is independent of the relaxation time across all regimes. Numerical experiments on applications in traffic flows and kinetic theory verify our theoretical analysis., Comment: arXiv admin note: text overlap with arXiv:2306.08742 by other authors

We introduce a hyperbolic closure for the Grad moment expansion of the Bhatnagar-Gross-Krook's (BGK) kinetic model using a neural network (NN) trained on BGK's moment data. This closure is motivated by the exact closure for the free streaming limit that we derived in our paper on closures in transport \cite{Huang2022-RTE1}. The exact closure relates the gradient of the highest moment to the gradient of four lower moments. As with our past work, the model presented here learns the gradient of the highest moment in terms of the coefficients of gradients for all lower ones. By necessity, this means that the resulting hyperbolic system is not conservative in the highest moment. For stability, the output layers of the NN are designed to enforce hyperbolicity and Galilean invariance. This ensures the model can be run outside of the training window of the NN. Unlike our previous work on radiation transport that dealt with linear models, the BGK model's nonlinearity demanded advanced training tools. These comprised an optimal learning rate discovery, one cycle training, batch normalization in each neural layer, and the use of the \texttt{AdamW} optimizer. To address the non-conservative structure of the hyperbolic model, we adopt the FORCE numerical method to achieve robust solutions. This results in a comprehensive computing model combining learned closures with methods for solving hyperbolic models. The proposed model can capture accurate moment solutions across a broad spectrum of Knudsen numbers. Our paper details the multi-scale model construction and is run on a range of test problems., Comment: 30 pages, 7 figures

This work is concerned with the uniform accuracy of implicit-explicit backward differentiation formulas for general linear hyperbolic relaxation systems satisfying the structural stability condition proposed previously by the third author. We prove the uniform stability and accuracy of a class of IMEX-BDF schemes discretized spatially by a Fourier spectral method. The result reveals that the accuracy of the fully discretized schemes is independent of the relaxation time in all regimes. It is verified by numerical experiments on several applications to traffic flows, rarefied gas dynamics and kinetic theory.

In this paper, we investigate the two-dimensional extension of a recently introduced set of shallow water models based on a regularized moment expansion of the incompressible Navier-Stokes equations \cite{kowalski2017moment,koellermeier2020analysis}. We show the rotational invariance of the proposed moment models with two different approaches. The first proof involves the split of the coefficient matrix into the conservative and non-conservative parts and proves the rotational invariance for each part, while the second one relies on the special block structure of the coefficient matrices. With the aid of rotational invariance, the analysis of the hyperbolicity for the moment model in 2D is reduced to the real diagonalizability of the coefficient matrix in 1D. Then we analyze the real diagonalizability by deriving the analytical form of the characteristic polynomial. We find that the moment model in 2D is hyperbolic in most cases and weakly hyperbolic in a degenerate edge case. With a simple modification to the coefficient matrices, we fix this weakly hyperbolicity and propose a new global hyperbolic model. Furthermore, we extend the model to include a more general class of closure relations than the original model and establish that this set of general closure relations retains both rotational invariance and hyperbolicity.

The flat band system is an ideal quantum platform to investigate the kaleidoscope created by the electron-electron correlation effects. The central ingredient of realizing a flat band is to find its compact localized states. In this work, we develop a systematic way to generate the compact localized states by designing destructive interference pattern from 1-dimensional chains. A variety of 2-dimensional new flat band systems are constructed with this method. Furthermore, we show that the method can be extended to generate the compact localized states in multi-orbital systems by carefully designing the block hopping scheme, as well as in quasicrystal and disorder systems.

In this paper, we develop bound-preserving discontinuous Galerkin (DG) methods for chemical reactive flows. There are several difficulties in constructing suitable numerical schemes. First of all, the density and internal energy are positive, and the mass fraction of each species is between 0 and 1. Secondly, due to the rapid reaction rate, the system may contain stiff sources, and the strong-stability-preserving explicit Runge-Kutta method may result in limited time step sizes. To obtain physically relevant numerical approximations, we apply the bound-preserving technique to the DG methods. For time discretization, we apply the modified Runge-Kutta/multi-step Patankar methods, which are explicit for the flux while implicit for the source. Such methods can handle stiff sources with relatively large time steps, preserve the positivity of the target variables, and keep the summation of the mass fractions to be 1. Finally, it is not straightforward to combine the bound-preserving DG methods and the Patankar time integrations. The positivity-preserving technique for DG method requires positive numerical approximations at the cell interfaces, while Patankar methods can keep the positivity of the pre-selected point-values of the target variables. To match the degree of freedom, we use $Q^k$ polynomials on rectangular meshes for problems in two space dimensions. To evolve in time, we first read the polynomials at the Gaussian points. Then suitable slope limiters can be applied to enforce the positivity of the solutions at those points, which can be preserved by the Patankar methods, leading to positive updated numerical cell averages. In addition, we use another slope limiter to get positive solutions used for the bound-preserving technique for the flux.

This paper reviews the adaptive sparse grid discontinuous Galerkin (aSG-DG) method for computing high dimensional partial differential equations (PDEs) and its software implementation. The C\texttt{++} software package called AdaM-DG, implementing the aSG-DG method, is available on Github at \url{https://github.com/JuntaoHuang/adaptive-multiresolution-DG}. The package is capable of treating a large class of high dimensional linear and nonlinear PDEs. We review the essential components of the algorithm and the functionality of the software, including the multiwavelets used, assembling of bilinear operators, fast matrix-vector product for data with hierarchical structures. We further demonstrate the performance of the package by reporting numerical error and CPU cost for several benchmark test, including linear transport equations, wave equations and Hamilton-Jacobi equations.

This paper considers the discontinuous Galerkin (DG) methods for solving the Vlasov-Maxwell (VM) system, a fundamental model for collisionless magnetized plasma. The DG methods provide accurate numerical description with conservation and stability properties. However, to resolve the high dimensional probability distribution function, the computational cost is the main bottleneck even for modern-day supercomputers. This work studies the applicability of a post-processing technique to the DG solution to enhance its accuracy and resolution for the VM system. In particular, we prove the superconvergence of order $(2k+\frac{1}{2})$ in the negative order norm for the probability distribution function and the electromagnetic fields when piecewise polynomial degree $k$ is used. Numerical tests including Landau damping, two-stream instability and streaming Weibel instabilities are considered showing the performance of the post-processor.

This work is concerned with coupling conditions for linear hyperbolic relaxation systems with multiple relaxation times. In the region with small relaxation time, an equilibrium system can be used for computational efficiency. Under the assumption that the relaxation system satisfies the structural stability condition and the interface is non-characteristic, we derive a coupling condition at the interface to couple the two systems in a domain decomposition setting. We prove the validity by the energy estimate and Laplace transform, which shows how the error of the domain decomposition method depends on the smaller relaxation time and the boundary layer effects. In addition, we propose a discontinuous Galerkin (DG) scheme for solving the interface problem with the derived coupling condition and prove the L2 stability. We validate our analysis on the linearized Carleman model and the linearized Grad's moment system and show the effectiveness of the DG scheme., Comment: 30 pages, 2 figures

In this paper, we perform stability analysis for a class of second and third order accurate strong-stability-preserving modified Patankar Runge-Kutta (SSPMPRK) schemes, which were introduced in [4,5] and can be used to solve convection equations with stiff source terms, such as reactive Euler equations, with guaranteed positivity under the standard CFL condition due to the convection terms only. The analysis allows us to identify the range of free parameters in these SSPMPRK schemes in order to ensure stability. Numerical experiments are provided to demonstrate the validity of the analysis., Comment: 23 pages, 10 figures

This is the third paper in a series in which we develop machine learning (ML) moment closure models for the radiative transfer equation (RTE). In our previous work \cite{huang2021gradient}, we proposed an approach to learn the gradient of the unclosed high order moment, which performs much better than learning the moment itself and the conventional $P_N$ closure. However, while the ML moment closure has better accuracy, it is not able to guarantee hyperbolicity and has issues with long time stability. In our second paper \cite{huang2021hyperbolic}, we identified a symmetrizer which leads to conditions that enforce that the gradient based ML closure is symmetrizable hyperbolic and stable over long time. The limitation of this approach is that in practice the highest moment can only be related to four, or fewer, lower moments. In this paper, we propose a new method to enforce the hyperbolicity of the ML closure model. Motivated by the observation that the coefficient matrix of the closure system is a lower Hessenberg matrix, we relate its eigenvalues to the roots of an associated polynomial. We design two new neural network architectures based on this relation. The ML closure model resulting from the first neural network is weakly hyperbolic and guarantees the physical characteristic speeds, i.e., the eigenvalues are bounded by the speed of light. The second model is strictly hyperbolic and does not guarantee the boundedness of the eigenvalues. Several benchmark tests including the Gaussian source problem and the two-material problem show the good accuracy, stability and generalizability of our hyperbolic ML closure model., Comment: arXiv admin note: text overlap with arXiv:2105.14410

This is the second paper in a series in which we develop machine learning (ML) moment closure models for the radiative transfer equation (RTE). In our previous work \cite{huang2021gradient}, we proposed an approach to directly learn the gradient of the unclosed high order moment, which performs much better than learning the moment itself and the conventional $P_N$ closure. However, the ML moment closure model in \cite{huang2021gradient} is not able to guarantee hyperbolicity and long time stability. We propose in this paper a method to enforce the global hyperbolicity of the ML closure model. The main idea is to seek a symmetrizer (a symmetric positive definite matrix) for the closure system, and derive constraints such that the system is globally symmetrizable hyperbolic. It is shown that the new ML closure system inherits the dissipativeness of the RTE and preserves the correct diffusion limit as the Knunsden number goes to zero. Several benchmark tests including the Gaussian source problem and the two-material problem show the good accuracy, long time stability and generalizability of our globally hyperbolic ML closure model.

In this paper, we take a data-driven approach and apply machine learning to the moment closure problem for radiative transfer equation in slab geometry. Instead of learning the unclosed high order moment, we propose to directly learn the gradient of the high order moment using neural networks. This new approach is consistent with the exact closure we derive for the free streaming limit and also provides a natural output normalization. A variety of benchmark tests, including the variable scattering problem, the Gaussian source problem with both periodic and reflecting boundaries, and the two-material problem, show both good accuracy and generalizability of our machine learning closure model.

In this paper, we propose a class of adaptive multiresolution (also called adaptive sparse grid) ultra-weak discontinuous Galerkin (UWDG) methods for solving some nonlinear dispersive wave equations including the Korteweg-de Vries (KdV) equation and its two dimensional generalization, the Zakharov-Kuznetsov (ZK) equation. The UWDG formulation, which relies on repeated integration by parts, was proposed for KdV equation in \cite{cheng2008discontinuous}. For the ZK equation which contains mixed derivative terms, we develop a new UWDG formulation. The $L^2$ stability and the optimal error estimate with a novel local projection are established for this new scheme on regular meshes. Adaptivity is achieved based on multiresolution and is particularly effective for capturing solitary wave structures. Various numerical examples are presented to demonstrate the accuracy and capability of our methods.

In this paper, we propose a data-driven method to discover multiscale chemical reactions governed by the law of mass action. First, we use a single matrix to represent the stoichiometric coefficients for both the reactants and products in a system without catalysis reactions. The negative entries in the matrix denote the stoichiometric coefficients for the reactants and the positive ones for the products. Second, we find that the conventional optimization methods usually get stuck in the local minima and could not find the true solution in learning the multiscale chemical reactions. To overcome this difficulty, we propose a partial-parameters-freezing (PPF) technique to progressively determine the network parameters by using the fact that the stoichiometric coefficients are integers. With such a technique, the dimension of the searching space is gradually reduced in the training process and the global mimina can be eventually obtained. Several numerical experiments including the classical Michaelis-Menten kinetics, the hydrogen oxidation reactions, and the simplified GRI-3.0 mechanism verify the good performance of our algorithm in learning the multiscale chemical reactions. The code is available at \url{https://github.com/JuntaoHuang/multiscale-chemical-reaction}.

In this work, we develop a method for learning interpretable, thermodynamically stable and Galilean invariant partial differential equations (PDEs) based on the Conservation-dissipation Formalism of irreversible thermodynamics. As governing equations for non-equilibrium flows in one dimension, the learned PDEs are parameterized by fully-connected neural networks and satisfy the conservation-dissipation principle automatically. In particular, they are hyperbolic balance laws and Galilean invariant. The training data are generated from a kinetic model with smooth initial data. Numerical results indicate that the learned PDEs can achieve good accuracy in a wide range of Knudsen numbers. Remarkably, the learned dynamics can give satisfactory results with randomly sampled discontinuous initial data and Sod's shock tube problem although it is trained only with smooth initial data.

This paper develops a high order adaptive scheme for solving nonlinear Schrodinger equations. The solutions to such equations often exhibit solitary wave and local structures, which makes adaptivity essential in improving the simulation efficiency. Our scheme uses the ultra-weak discontinuous Galerkin (DG) formulation and belongs to the framework of adaptive multiresolution schemes. Various numerical experiments are presented to demonstrate the excellent capability of capturing the soliton waves and the blow-up phenomenon.

We are interested in numerically solving the Hamilton-Jacobi (HJ) equations, which arise in optimal control and many other applications. Oftentimes, such equations are posed in high dimensions, and this poses great numerical challenges. This work proposes a class of adaptive sparse grid (also called adaptive multiresolution) local discontinuous Galerkin (DG) methods for solving Hamilton-Jacobi equations in high dimensions. By using the sparse grid techniques, we can treat moderately high dimensional cases. Adaptivity is incorporated to capture kinks and other local structures of the solutions. Two classes of multiwavelets are used to achieve multiresolution, which are the orthonormal Alpert's multiwavelets and the interpolatory multiwavelets. Numerical tests in up to four dimensions are provided to validate the performance of the method., Comment: arXiv admin note: text overlap with arXiv:2004.08525

In this paper, we propose a class of adaptive multiresolution (also called adaptive sparse grid) discontinuous Galerkin (DG) methods for simulating scalar wave equations in second order form in space. The two key ingredients of the schemes include an interior penalty DG formulation in the adaptive function space and two classes of multiwavelets for achieving multiresolution. In particular, the orthonormal Alpert's multiwavelets are used to express the DG solution in terms of a hierarchical structure, and the interpolatory multiwavelets are further introduced to enhance computational efficiency in the presence of variable wave speed or nonlinear source. Some theoretical results on stability and accuracy of the proposed method are presented. Benchmark numerical tests in 2D and 3D are provided to validate the performance of the method.

In \cite{ZH2019}, we developed a boundary treatment method for implicit-explicit (IMEX) Runge-Kutta (RK) methods for solving hyperbolic systems with source terms. Since IMEX RK methods include explicit ones as special cases, this boundary treatment method naturally applies to explicit methods as well. In this paper, we examine this boundary treatment method for the case of explicit RK schemes of arbitrary order applied to hyperbolic conservation laws. We show that the method not only preserves the accuracy of explicit RK schemes but also possesses good stability. This compares favourably to the inverse Lax-Wendroff method in \cite{TS2010,TWSN2012} where analysis and numerical experiments have previously verified the presence of order reduction \cite{TS2010,TWSN2012}. In addition, we demonstrate that our method performs well for strong-stability-preserving (SSP) RK schemes involving negative coefficients and downwind spatial discretizations. It is numerically shown that when boundary conditions are present and the proposed boundary treatment is used, that SSP RK schemes with negative coefficients still allow for larger time steps than schemes with all non-negative coefficients. In this regard, our boundary treatment method is an effective supplement to SSP RK schemes with/without negative coefficients for initial-boundary value problems for hyperbolic conservation laws., Comment: To appear in Journal of Computational Physics. arXiv admin note: text overlap with arXiv:1908.01027

In this paper, we develop a high order finite difference boundary treatment method for the implicit-explicit (IMEX) Runge-Kutta (RK) schemes solving hyperbolic systems with possibly stiff source terms on a Cartesian mesh. The main challenge is how to obtain the solutions at ghost points resulting from the wide stencil of the interior high order scheme. We address this problem by combining the idea of using the RK schemes at the boundary and an inverse Lax-Wendroff procedure. The former preserves the accuracy of the RK schemes and the latter guarantees the stability. Our method is different from the widely used approach for the explicit RK schemes by imposing boundary conditions at intermediate stages, which could not be derived for the IMEX schemes. In addition, the intermediate boundary conditions are only available for explicit RK schemes up to third order while our method applies to arbitrary order IMEX and explicit RK schemes. Moreover, the present boundary treatment method may be adapted to IMEX RK schemes solving many other partial differential equations. For a specific third-order IMEX scheme, we demonstrate the good stability and third-order accuracy of our boundary treatment through both 1D examples and 2D reactive Euler equations.

In this paper, we develop an adaptive multiresolution discontinuous Galerkin (DG) scheme for scalar hyperbolic conservation laws in multidimensions. Compared with previous work for linear hyperbolic equations \cite{guo2016transport, guo2017adaptive}, a class of interpolatory multiwavelets are applied to efficiently compute the nonlinear integrals over elements and edges in DG schemes. The resulting algorithm, therefore can achieve similar computational complexity as the sparse grid DG method for smooth solutions. Theoretical and numerical studies are performed taking into consideration of accuracy and stability with regard to the choice of the interpolatory multiwavelets. Artificial viscosity is added to capture the shock and only acts on the leaf elements taking advantages of the multiresolution representation. Adaptivity is realized by auto error thresholding based on hierarchical surplus. Accuracy and robustness are demonstrated by several numerical tests.