Search

Time-dependent kinetic models are ubiquitous in computational science and engineering. The underlying integro-differential equations in these models are high-dimensional, comprised of a six--dimensional phase space, making simulations of such phenomena extremely expensive. In this article we demonstrate that in many situations, the solution to kinetics problems lives on a low dimensional manifold that can be described by a low-rank matrix or tensor approximation. We then review the recent development of so-called low-rank methods that evolve the solution on this manifold. The two classes of methods we review are the dynamical low-rank (DLR) method, which derives differential equations for the low-rank factors, and a Step-and-Truncate (SAT) approach, which projects the solution onto the low-rank representation after each time step. Thorough discussions of time integrators, tensor decompositions, and method properties such as structure preservation and computational efficiency are included. We further show examples of low-rank methods as applied to particle transport and plasma dynamics.

High-order semi-Lagrangian methods for kinetic equations have been under rapid development in the past few decades. In this work, we propose a semi-Lagrangian adaptive rank (SLAR) integrator in the finite difference framework for linear advection and nonlinear Vlasov-Poisson systems without dimensional splitting. The proposed method leverages the semi-Lagrangian approach to allow for significantly larger time steps while also exploiting the low-rank structure of the solution. This is achieved through cross approximation of matrices, also referred to as CUR or pseudo-skeleton approximation, where representative columns and rows are selected using specific strategies. To maintain numerical stability and ensure local mass conservation, we apply singular value truncation and a mass-conservative projection following the cross approximation of the updated solution. The computational complexity of our method scales linearly with the mesh size $N$ per dimension, compared to the $\mathcal{O}(N^2)$ complexity of traditional full-rank methods per time step. The algorithm is extended to handle nonlinear Vlasov-Poisson systems using a Runge-Kutta exponential integrator. Moreover, we evolve the macroscopic conservation laws for charge densities implicitly, enabling the use of large time steps that align with the semi-Lagrangian solver. We also perform a mass-conservative correction to ensure that the adaptive rank solution preserves macroscopic charge density conservation. To validate the efficiency and effectiveness of our method, we conduct a series of benchmark tests on both linear advection and nonlinear Vlasov-Poisson systems. The propose algorithm will have the potential in overcoming the curse of dimensionality for beyond 2D high dimensional problems, which is the subject of our future work., Comment: 25 pages, 15 figures

We consider the adaptive-rank integration of general time-dependent advection-diffusion partial differential equations (PDEs) with spatially variable coefficients. We employ a standard finite-difference method for spatial discretization coupled with diagonally implicit Runge-Kutta schemes for temporal discretization. The fully discretized scheme can be written as a generalized Sylvester equation, which we solve with an adaptive-rank algorithm structured around three key strategies: (i) constructing dimension-wise subspaces based on an extended Krylov strategy, (ii) developing an effective Averaged-Coefficient Sylvester (ACS) preconditioner to invert the reduced system for the coefficient matrix efficiently, and (iii) efficiently computing the residual of the equation without explicitly reverting to the full-rank form. The proposed approach features a computational complexity of O(Nr2 + r3), where r represents the rank during the Krylov iteration and N is the resolution in one dimension, commensurate with the constant-coefficient case [El Kahza et al, J. Comput. Phys., 518 (2024)]. We present numerical examples that illustrate the computational efficacy and complexity of our algorithm.

The tensor-train (TT) format is a data-sparse tensor representation commonly used in high dimensional function approximations arising from computational and data sciences. Various sequential and parallel TT decomposition algorithms have been proposed for different tensor inputs and assumptions. In this paper, we propose subtensor parallel adaptive TT cross, which partitions a tensor onto distributed memory machines with multidimensional process grids, and constructs an TT approximation iteratively with tensor elements. We derive two iterative formulations for pivot selection and TT core construction under the distributed memory setting, conduct communication and scaling analysis of the algorithm, and illustrate its performance with multiple test experiments. These include up to 6D Hilbert tensors and tensors constructed from Maxwellian distribution functions that arise in kinetic theory. Our results demonstrate significant accuracy with greatly reduced storage requirements via the TT cross approximation. Furthermore, we demonstrate good to optimal strong and weak scaling performance for the proposed parallel algorithm.

In this paper, we present a new adaptive rank approximation technique for computing solutions to the high-dimensional linear kinetic transport equation. The approach we propose is based on a macro-micro decomposition of the kinetic model in which the angular domain is discretized with a tensor product quadrature rule under the discrete ordinates method. To address the challenges associated with the curse of dimensionality, the proposed low-rank method is cast in the framework of the hierarchical Tucker decomposition. The adaptive rank integrators we propose are built upon high-order discretizations for both time and space. In particular, this work considers implicit-explicit discretizations for time and finite-difference weighted-essentially non-oscillatory discretizations for space. The high-order singular value decomposition is used to perform low-rank truncation of the high-dimensional time-dependent distribution function. The methods are applied to several benchmark problems, where we compare the solution quality and measure compression achieved by the adaptive rank methods against their corresponding full-grid methods. We also demonstrate the benefits of high-order discretizations in the proposed low-rank framework., Comment: 26 pages, 15 figures, 2 tables, 48 references

In this paper, we develop high-order, conservative, non-splitting Eulerian-Lagrangian (EL) Runge-Kutta (RK) finite volume (FV) weighted essentially non-oscillatory (WENO) schemes for convection-diffusion equations. The proposed EL-RK-FV-WENO scheme defines modified characteristic lines and evolves the solution along them, significantly relaxing the time-step constraint for the convection term. The main algorithm design challenge arises from the complexity of constructing accurate and robust reconstructions on dynamically varying Lagrangian meshes. This reconstruction process is needed for flux evaluations on time-dependent upstream quadrilaterals and time integrations along moving characteristics. To address this, we propose a strategy that utilizes a WENO reconstruction on a fixed Eulerian mesh for spatial reconstruction, and updates intermediate solutions on the Eulerian background mesh for implicit-explicit RK temporal integration. This strategy leverages efficient reconstruction and remapping algorithms to manage the complexities of polynomial reconstructions on time-dependent quadrilaterals, while ensuring local mass conservation. The proposed scheme ensures mass conservation due to the flux-form semi-discretization and the mass-conservative reconstruction on both background and upstream cells. Extensive numerical tests have been performed to verify the effectiveness of the proposed scheme.

We present a class of high-order Eulerian-Lagrangian Runge-Kutta finite volume methods that can numerically solve Burgers' equation with shock formations, which could be extended to general scalar conservation laws. Eulerian-Lagrangian (EL) and semi-Lagrangian (SL) methods have recently seen increased development and have become a staple for allowing large time-stepping sizes. Yet, maintaining relatively large time-stepping sizes post shock formation remains quite challenging. Our proposed scheme integrates the partial differential equation on a space-time region partitioned by linear approximations to the characteristics determined by the Rankine-Hugoniot jump condition. We trace the characteristics forward in time and present a merging procedure for the mesh cells to handle intersecting characteristics due to shocks. Following this partitioning, we write the equation in a time-differential form and evolve with Runge-Kutta methods in a method-of-lines fashion. High-resolution methods such as ENO and WENO-AO schemes are used for spatial reconstruction. Extension to higher dimensions is done via dimensional splitting. Numerical experiments demonstrate our scheme's high-order accuracy and ability to sharply capture post-shock solutions with large time-stepping sizes., Comment: 29 pages

We propose a high order adaptive-rank implicit integrators for stiff time-dependent PDEs, leveraging extended Krylov subspaces to efficiently and adaptively populate low-rank solution bases. This allows for the accurate representation of solutions with significantly reduced computational costs. We further introduce an efficient mechanism for residual evaluation and an adaptive rank-seeking strategy that optimizes low-rank settings based on a comparison between the residual size and the local truncation errors of the time-stepping discretization. We demonstrate our approach with the challenging Lenard-Bernstein Fokker-Planck (LBFP) nonlinear equation, which describes collisional processes in a fully ionized plasma. The preservation of {the equilibrium state} is achieved through the Chang-Cooper discretization, and strict conservation of mass, momentum and energy via a Locally Macroscopic Conservative (LoMaC) procedure. The development of implicit adaptive-rank integrators, demonstrated here up to third-order temporal accuracy via diagonally implicit Runge-Kutta schemes, showcases superior performance in terms of accuracy, computational efficiency, equilibrium preservation, and conservation of macroscopic moments. This study offers a starting point for developing scalable, efficient, and accurate methods for high-dimensional time-dependent problems., Comment: 29 pages, 8 figures

This paper introduces a novel computational approach termed the Reduced Augmentation Implicit Low-rank (RAIL) method by investigating two predominant research directions in low-rank solutions to time-dependent partial differential equations (PDEs): dynamical low-rank (DLR), and step and truncation (SAT) tensor methods. The RAIL method, along with the development of the SAT approach, is designed to enhance the efficiency of traditional full-rank implicit solvers from method-of-lines discretizations of time-dependent PDEs, while maintaining accuracy and stability. We consider spectral methods for spatial discretization, and diagonally implicit Runge-Kutta (DIRK) and implicit-explicit (IMEX) RK methods for time discretization. The efficiency gain is achieved by investigating low-rank structures within solutions at each RK stage using a singular value decomposition (SVD). In particular, we develop a reduced augmentation procedure to predict the basis functions to construct projection subspaces. This procedure balances algorithm accuracy and efficiency by incorporating as many bases as possible from previous RK stages and predictions, and by optimizing the basis representation through SVD truncation. As such, one can form implicit schemes for updating basis functions in a dimension-by-dimension manner, similar in spirit to the K-L step in the DLR framework. We also apply a globally mass conservative post-processing step at the end of each RK stage. We validate the RAIL method through numerical simulations of advection-diffusion problems and a Fokker-Planck model, showcasing its ability to efficiently handle time-dependent PDEs while maintaining global mass conservation. Our approach generalizes and bridges the DLR and SAT approaches, offering a comprehensive framework for efficiently and accurately solving time-dependent PDEs with implicit treatment.

In this paper, we construct a novel Eulerian-Lagrangian finite volume (ELFV) method for nonlinear scalar hyperbolic equations in one space dimension. It is well known that the exact solutions to such problems may contain shocks though the initial conditions are smooth, and direct numerical methods may suffer from restricted time step sizes. To relieve the restriction, we propose an ELFV method, where the space-time domain was separated by the partition lines originated from the cell interfaces whose slopes are obtained following the Rakine-Hugoniot junmp condition. Unfortunately, to avoid the intersection of the partition lines, the time step sizes are still limited. To fix this gap, we detect effective troubled cells (ETCs) and carefully design the influence region of each ETC, within which the partitioned space-time regions are merged together to form a new one. Then with the new partition of the space-time domain, we theoretically prove that the proposed first-order scheme with Euler forward time discretization is total-variation-diminishing and maximum-principle-preserving with {at least twice} larger time step constraints than the classical first order Eulerian method for Burgers' equation. Numerical experiments verify the optimality of the designed time step sizes., Comment: Submitted to Mathematics of Computation

In this paper, we propose a novel Local Macroscopic Conservative (LoMaC) low rank tensor method with discontinuous Galerkin (DG) discretization for the physical and phase spaces for simulating the Vlasov-Poisson (VP) system. The LoMaC property refers to the exact local conservation of macroscopic mass, momentum and energy at the discrete level. The recently developed LoMaC low rank tensor algorithm (arXiv:2207.00518) simultaneously evolves the macroscopic conservation laws of mass, momentum and energy using the kinetic flux vector splitting; then the LoMaC property is realized by projecting the low rank kinetic solution onto a subspace that shares the same macroscopic observables. This paper is a generalization of our previous work, but with DG discretization to take advantage of its compactness and flexibility in handling boundary conditions and its superior accuracy in the long term. The algorithm is developed in a similar fashion as that for a finite difference scheme, by observing that the DG method can be viewed equivalently in a nodal fashion. With the nodal DG method, assuming a tensorized computational grid, one will be able to (1) derive differentiation matrices for different nodal points based on a DG upwind discretization of transport terms, and (2) define a weighted inner product space based on the nodal DG grid points. The algorithm can be extended to the high dimensional problems by hierarchical Tucker decomposition of solution tensors and a corresponding conservative projection algorithm. In a similar spirit, the algorithm can be extended to DG methods on nodal points of an unstructured mesh, or to other types of discretization, e.g. the spectral method in velocity direction. Extensive numerical results are performed to showcase the efficacy of the method., Comment: arXiv admin note: text overlap with arXiv:2207.00518

We present fourth-order conservative non-splitting semi-Lagrangian (SL) Hermite essentially non-oscillatory (HWENO) schemes for linear transport equations with applications for nonlinear problems including the Vlasov-Poisson system, the guiding center Vlasov model, and the incompressible Euler equations in the vorticity-stream function formulation. The proposed SL HWENO schemes combine a weak formulation of the characteristic Galerkin method with two newly constructed HWENO reconstruction methods. Fourth-order accuracy is accomplished in both space and time under a non-splitting setting. Mass conservation naturally holds due to the weak formulation of the characteristic Galerkin method and the design of the HWENO reconstructions. We apply a positive-preserving limiter to maintain the positivity of numerical solutions when needed. Although the proposed SL framework allows us to take large time steps for improving computational efficiency, it also brings challenges to the spatial reconstruction technique; we construct two kind of novel HWENO reconstructions to fit the need for the proposed SL framework. Abundant benchmark tests are performed to verify the effectiveness of the proposed SL HWENO schemes.

We propose an Eulerian-Lagrangian (EL) Runge-Kutta (RK) discontinuous Galerkin (DG) method for wave equations. The method is designed based on the ELDG method for transport problems [J. Comput. Phy. 446: 110632, 2021.], which tracks solution along approximations to characteristics in the DG framework, allowing extra large time stepping sizes with stability. The wave equation can be written as a first order hyperbolic system. Considering each characteristic family, a straightforward application of ELDG will be to transform to the characteristic variables, evolve them on associated characteristic related space-time regions, and transform them back to the original variables. However, the mass conservation could not be guaranteed in a general setting. In this paper, we formulate a mass conservative semi-discrete ELDG method by decomposing each variable into two parts, each of them associated with a different characteristic family. As a result, four different quantities are evolved in EL fashion and recombined to update the solution. The fully discrete scheme is formulated by using method-of-lines RK methods, with intermediate RK solutions updated on the background mesh. Numerical results on 1D and 2D wave equations are presented to demonstrate the performance of the proposed ELDG method. These include the high order spatial and temporal accuracy, stability with extra large time stepping size, and mass conservative property., Comment: 31 pages

In this paper, we propose a novel Local Macroscopic Conservative (LoMaC) low rank tensor method for simulating the Vlasov-Poisson (VP) system. The LoMaC property refers to the exact local conservation of macroscopic mass, momentum and energy at the discrete level. This is a follow-up work of our previous development of a conservative low rank tensor approach for Vlasov dynamics (arXiv:2201.10397). In that work, we applied a low rank tensor method with a conservative singular value decomposition (SVD) to the high dimensional VP system to mitigate the curse of dimensionality, while maintaining the local conservation of mass and momentum. However, energy conservation is not guaranteed, which is a critical property to avoid unphysical plasma self-heating or cooling. The new ingredient in the LoMaC low rank tensor algorithm is that we simultaneously evolve the macroscopic conservation laws of mass, momentum and energy using a flux-difference form with kinetic flux vector splitting; then the LoMaC property is realized by projecting the low rank kinetic solution onto a subspace that shares the same macroscopic observables by a conservative orthogonal projection. The algorithm is extended to the high dimensional problems by hierarchical Tuck decomposition of solution tensors and a corresponding conservative projection algorithm. Extensive numerical tests on the VP system are showcased for the algorithm's efficacy., Comment: arXiv admin note: substantial text overlap with arXiv:2201.10397

In this paper, we propose a conservative low rank tensor method to approximate nonlinear Vlasov solutions. The low rank approach is based on our earlier work (arxiv: 2106.08834). It takes advantage of the fact that the differential operators in the Vlasov equation are tensor friendly, based on which we propose to dynamically and adaptively build up low rank solution basis by adding new basis functions from discretization of the differential equation, and removing basis from a singular value decomposition (SVD)-type truncation procedure. For the discretization, we adopt a high order finite difference spatial discretization together with a second order strong stability preserving multi-step time discretization. While the SVD truncation will remove the redundancy in representing the high dimensional Vlasov solution, it will destroy the conservation properties of the associated full conservative scheme. In this paper, we develop a conservative truncation procedure with conservation of mass, momentum and kinetic energy densities. The conservative truncation is achieved by an orthogonal projection onto a subspace spanned by $1$, $v$ and $v^2$ in the velocity space associated with a weighted inner product. Then the algorithm performs a weighted SVD truncation of the remainder, which involves a scaling, followed by the standard SVD truncation and rescaling back. The algorithm is further developed in high dimensions with hierarchical Tucker tensor decomposition of high dimensional Vlasov solutions, overcoming the curse of dimensionality. An extensive set of nonlinear Vlasov examples are performed to show the effectiveness and conservation property of proposed conservative low rank approach. Comparison is performed against the non-conservative low rank tensor approach on conservation history of mass, momentum and energy.

We propose a low-rank tensor approach to approximate linear transport and nonlinear Vlasov solutions and their associated flow maps. The approach takes advantage of the fact that the differential operators in the Vlasov equation is tensor friendly, based on which we propose a novel way to dynamically and adaptively build up low-rank solution basis by adding new basis functions from discretization of the PDE, and removing basis from an SVD-type truncation procedure. For the discretization, we adopt a high order finite difference spatial discretization and a second order strong stability preserving multi-step time discretization. We apply the same procedure to evolve the dynamics of the flow map in a low-rank fashion, which proves to be advantageous when the flow map enjoys the low rank structure, while the solution suffers from high rank or displays filamentation structures. Hierarchical Tucker decomposition is adopted for high dimensional problems. An extensive set of linear and nonlinear Vlasov test examples are performed to show the high order spatial and temporal convergence of the algorithm with mesh refinement up to SVD-type truncation, the significant computational savings of the proposed low-rank approach especially for high dimensional problems, the improved performance of the flow map approach for solutions with filamentations.

In this paper, we propose a new high order semi-implicit scheme for the all Mach full Euler equations of gas dynamics. Material waves are treated explicitly, while acoustic waves are treated implicitly, thus avoiding severe CFL restrictions for low Mach flows. High order accuracy in time is obtained by semi-implicit temporal integrator based on the IMEX Runge-Kutta (IMEX-RK) framework. High order in space is achieved by finite difference WENO schemes with characteristic-wise reconstructions adapted to the semi-implicit IMEX-RK time discretization. Type A IMEX schemes are constructed to handle not well-prepared initial conditions. Besides, these schemes are proven to be asymptotic preserving and asymptotically accurate as the Mach number vanishes for well-prepared initial conditions. Divergence-free property of the time-discrete schemes is proved. The proposed scheme can also well capture discontinuous solutions in the compressible regime, especially for two dimensional Riemann problems. Numerical tests in one and two space dimensions will illustrate the effectiveness of the proposed schemes.

In this paper, we develop a family of third order asymptotic-preserving (AP) and asymptotically accurate (AA) diagonally implicit Runge-Kutta (DIRK) time discretization methods for the stiff hyperbolic relaxation systems and kinetic Bhatnagar-Gross-Krook (BGK) model in the semi-Lagrangian (SL) setting. The methods are constructed based on an accuracy analysis of the SL scheme for stiff hyperbolic relaxation systems and kinetic BGK model in the limiting fluid regime when the Knudsen number approaches $0$. An extra order condition for the asymptotic third order accuracy in the limiting regime is derived. Linear Von Neumann stability analysis of the proposed third order DIRK methods are performed to a simplified two-velocity linear kinetic model. Extensive numerical tests are presented to demonstrate the AA, AP and stability properties of our proposed schemes.

In this paper, we propose an efficient, high order accurate and asymptotic-preserving (AP) semi-Lagrangian (SL) method for the BGK model with constant or spatially dependent Knudsen number. The spatial discretization is performed by a mass conservative nodal discontinuous Galerkin (NDG) method, while the temporal discretization of the stiff relaxation term is realized by stiffly accurate diagonally implicit Runge-Kutta (DIRK) methods along characteristics. Extra order conditions are enforced for asymptotic accuracy (AA) property of DIRK methods when they are coupled with a semi-Lagrangian algorithm in solving the BGK model. A local maximum principle preserving (LMPP) limiter is added to control numerical oscillations in the transport step. Thanks to the SL and implicit nature of time discretization, the time stepping constraint is relaxed and it is much larger than that from an Eulerian framework with explicit treatment of the source term. Extensive numerical tests are presented to verify the high order AA, efficiency and shock capturing properties of the proposed schemes.

We propose a generalized Eulerian-Lagrangian (GEL) discontinuous Galerkin (DG) method. The method is a generalization of the Eulerian-Lagrangian (EL) DG method for transport problems proposed in [arXiv preprint arXiv: 2002.02930 (2020)], which tracks solution along approximations to characteristics in the DG framework, allowing extra large time stepping size with stability. The newly proposed GEL DG method in this paper is motivated for solving linear hyperbolic systems with variable coefficients, where the velocity field for adjoint problems of the test functions is frozen to constant. In this paper, in a simplified scalar setting, we propose the GEL DG methodology by freezing the velocity field of adjoint problems, and by formulating the semi-discrete scheme over the space-time region partitioned by linear lines approximating characteristics. The fully-discrete schemes are obtained by method-of-lines Runge-Kutta methods. We further design flux limiters for the schemes to satisfy the discrete geometric conservation law (DGCL) and maximum principle preserving (MPP) properties. Numerical results on 1D and 2D linear transport problems are presented to demonstrate great properties of the GEL DG method. These include the high order spatial and temporal accuracy, stability with extra large time stepping size, and satisfaction of DGCL and MPP properties.

In our recent work [22], a family of high order asymptotic preserving (AP) methods, termed as IMEX-LDG methods, are designed to solve some linear kinetic transport equations, including the one-group transport equation in slab geometry and the telegraph equation, in a diffusive scaling. As the Knudsen number $\varepsilon$ goes to zero, the limiting schemes are implicit discretizations to the limiting diffusive equation. Both Fourier analysis and numerical experiments imply the methods are unconditionally stable in the diffusive regime when $\varepsilon\ll1$. In this paper, we develop an energy approach to establish the numerical stability of the IMEX1-LDG method, the sub-family of the methods that is first order accurate in time and arbitrary order in space, for the model with general material properties. Our analysis is the first to simultaneously confirm unconditional stability when $\varepsilon\ll1$ and the uniform stability property with respect to $\varepsilon$. To capture the unconditional stability, a novel discrete energy is introduced by better exploring the contribution of the scattering term in different regimes. A general form of the weight function, introduced to obtain the unconditional stability for $\varepsilon\ll1$, is also for the first time considered in such stability analysis. Based on the uniform stability, a rigorous asymptotic analysis is then carried out to show the AP property., Comment: 24 pages, 1 figure

We propose a new Eulerian-Lagrangian (EL) discontinuous Galerkin (DG) method. The method is designed as a generalization of the semi-Lagrangian (SL) DG method for linear advection problems proposed in [J. Sci. Comput. 73: 514-542, 2017], which is formulated based on an adjoint problem and tracing upstream cells by tracking characteristics curves highly accurately. In the SLDG method, depending on the velocity field, upstream cells could be of arbitrary shape. Thus, a more sophisticated approximation to sides of the upstream cells is required to get high order approximation. For example, quadratic-curved (QC) quadrilaterals were proposed to approximate upstream cells for a third-order spatial accuracy in a swirling deformation example. In this paper, for linear advection problems, we propose a more general formulation, named the ELDG method. The scheme is formulated based on a {\em modified} adjoint problem for which the upstream cells are always quadrilaterals, which avoids the need to use QC quadrilaterals in the SLDG algorithm. The newly proposed ELDG method can be viewed as a new general framework, in which both the classical Eulerian Runge-Kutta DG formulation and the SL DG formulation can fit in. Numerical results on linear transport problems, as well as the nonlinear Vlasov and incompressible Euler dynamics using the exponential RK time integrators, are presented to demonstrate the effectiveness of the ELDG method.

In this paper, we propose a semi-Lagrangian discontinuous Galerkin method coupled with Runge-Kutta exponential integrators (SLDG-RKEI) for nonlinear Vlasov dynamics. The commutator-free Runge-Kutta (RK) exponential integrators (EI) were proposed by Celledoni, et al. (FGCS, 2003). In the nonlinear transport setting, the RKEI can be used to decompose the evolution of the nonlinear transport into a composition of a sequence of linearized dynamics. The resulting linearized transport equations can be solved by the semi-Lagrangian (SL) discontinuous Galerkin (DG) method proposed in Cai, et al. (JSC, 2017). The proposed method can achieve high order spatial accuracy via the SLDG framework, and high order temporal accuracy via the RK EI. Due to the SL nature, the proposed SLDG-RKEI method is not subject to the CFL condition, thus they have the potential in using larger time-stepping sizes than those in the Eulerian approach. Inheriting advantages from the SLDG method, the proposed SLDG-RKEI schemes are mass conservative, positivity-preserving, have no dimensional splitting error, perform well in resolving complex solution structures, and can be evolved with adaptive time-stepping sizes. We show the performance of the SLDG-RKEI algorithm by classical test problems for the nonlinear Vlasov-Poisson system, as well as the Guiding center Vlasov model. Though that it is not our focus of this paper to explore the SLDG-RKEI scheme for nonlinear hyperbolic conservation laws that develop shocks, we show some preliminary results on schemes' performance on the Burgers' equation.

In this paper, we propose an efficient high order semi-Lagrangian (SL) discontinuous Galerkin (DG) method for solving linear convection-diffusion equations. The method generalizes our previous work on developing the SLDG method for transport equations (J. Sci. Comput. 73: 514-542, 2017), making it capable of handling additional diffusion and source terms. Within the DG framework, the solution is evolved along the characteristics; while the diffusion term is discretized by the local DG (LDG) method and integrated along characteristics by implicit Runge-Kutta methods together with source terms. The proposed method is named the `SLDG-LDG' method and enjoys many attractive features of the DG and SL methods. These include the uniformly high order accuracy (e.g. third order) in space and in time, compact, mass conservative, and stability under large time stepping size. An $L^2$ stability analysis is provided when the method is coupled with the first order backward Euler discretization. Effectiveness of the method are demonstrated by a group of numerical tests in one and two dimensions., Comment: 34 pages; 16 figures; Journal of Computational Physics, accepted

Transport problems arise across diverse fields of science and engineering. Semi-Lagrangian (SL) discontinuous Galerkin (DG) methods are a class of high order deterministic transport solvers that enjoy advantages of both SL approach and DG spatial discretization. In this paper, we review existing SLDG methods to date and compare numerical their performances. In particular, we make a comparison between the splitting and non-splitting SLDG methods for multi-dimensional transport simulations. Through extensive numerical results, we offer a practical guide for choosing optimal SLDG solvers for linear and nonlinear transport simulations.

In this paper, we develop a class of high order conservative semi-Lagrangian (SL) discontinuous Galerkin (DG) methods for solving multi-dimensional linear transport equations. The methods rely on a characteristic Galerkin weak formulation, leading to $L^2$ stable discretizations for linear problems. Unlike many existing SL methods, the high order accuracy and mass conservation of the proposed methods are realized in a non-splitting manner. Thus, the detrimental splitting error, which is known to significantly contaminate long term transport simulations, will be not incurred. One key ingredient in the scheme formulation, borrowed from CSLAM [Lauritzen, Nair, and Ullrich, 2010], is the use of Green's theorem which allows us to convert volume integrals into a set of line integrals. The resulting line integrals are much easier to approximate with high order accuracy, hence facilitating the implementation. Another novel ingredient is the construction of quadratic curves in approximating sides of upstream cell, leading to quadratic-curved quadrilateral upstream cells. Formal third order accuracy is obtained by such a construction. The desired positivity-preserving property is further attained by incorporating a high order bound-preserving filter. To assess the performance of the proposed methods, we test and compare the numerical schemes with a variety of configurations for solving several benchmark transport problems with both smooth and nonsmooth solutions. The efficiency and efficacy are numerically verified.

In this paper, we propose a mass conservative semi-Lagrangian finite difference scheme for multi-dimensional problems without dimensional splitting. The semi-Lagrangian scheme, based on tracing characteristics backward in time from grid points, does not necessarily conserve the total mass. To ensure mass conservation, we propose a conservative correction procedure based on a flux difference form. Such procedure guarantees local mass conservation, while introducing time step constraints for stability. We theoretically investigate such stability constraints from an ODE point of view by assuming exact evaluation of spatial differential operators and from the Fourier analysis for linear PDEs. The scheme is tested by classical two dimensional linear passive-transport problems, such as linear advection, rotation and swirling deformation. The scheme is applied to solve the nonlinear Vlasov-Poisson system using a a high order tracing mechanism proposed in [Qiu and Russo, 2016]. Such high order characteristics tracing scheme is generalized to the nonlinear guiding center Vlasov model and incompressible Euler system. The effectiveness of the proposed conservative semi-Lagrangian scheme is demonstrated numerically by our extensive numerical tests.

In this paper, we consider a finite difference grid-based semi-Lagrangian approach in solving the Vlasov-Poisson (VP) system. Many of existing methods are based on dimensional splitting, which decouples the problem into solving linear advection problems, see {\em Cheng and Knorr, Journal of Computational Physics, 22(1976)}. However, such splitting is subject to the splitting error. If we consider multi-dimensional problems without splitting, difficulty arises in tracing characteristics with high order accuracy. Specifically, the evolution of characteristics is subject to the electric field which is determined globally from the distribution of particle densities via the Poisson's equation. In this paper, we propose a novel strategy of tracing characteristics high order in time via a two-stage multi-derivative prediction-correction approach and by using moment equations of the VP system. With the foot of characteristics being accurately located, we proposed to use weighted essentially non-oscillatory (WENO) interpolation to recover function values between grid points, therefore to update solutions at the next time level. The proposed algorithm does not have time step restriction as Eulerian approach and enjoys high order spatial and temporal accuracy. However, such finite difference algorithm does not enjoy mass conservation; we discuss one possible way of resolving such issue and its potential challenge in numerical stability. The performance of the proposed schemes are numerically demonstrated via classical test problems such as Landau damping and two stream instabilities.

In this paper, we present a high order conservative semi-Lagrangian (SL) Hermite weighted essentially non-oscillatory (HWENO) method for the Vlasov equation based on dimensional splitting [Cheng and Knorr, Journal of Computational Physics, 22(1976)]. The major advantage of HWENO reconstruction, compared with the original WENO reconstruction, is compact. For the split one-dimensional equation, to ensure local mass conservation, we propose a high order SL HWENO scheme in a conservative flux-difference form, following the work in [J.-M. Qiu and A. Christlieb, Journal of Computational Physics, v229(2010)]. Besides performing dimensional splitting for the original 2D problem, we design a proper splitting for equations of derivatives to ensure local mass conservation of the proposed HWENO scheme. The proposed fifth order SL HWENO scheme with the Eulerian CFL condition has been tested to work well in capturing filamentation structures without introducing oscillations. We introduce WENO limiters to control oscillations when the time stepping size is larger than the Eulerian CFL restriction. We perform classical numerical tests on rigid body rotation problem, and demonstrate the performance of our scheme via the Landau damping and two-stream instabilities when solving the Vlasov-Poisson system.

We illustrate that numerical solutions of high order finite volume Hermite weighted essentially non-oscillatory (HWENO) scheme for some nonconvex conservation laws perform poorly or converge to the entropy solution in a slow speed. The modified finite volume HWENO schemes based either on first order monotone schemes or a second order entropic projection following the work of Qiu and Shu [SIAM J. Sci. Comput., 31 (2008), 584-607] are proposed and compared for solving one-dimensional scalar problems. We extend the modified finite volume HWENO based on first order monotone schemes for one-dimensional systems and two-dimensional scalar conservation laws. Numerical tests for several representative examples will be reported.

In this paper, we develop high-order asymptotic preserving (AP) schemes for the BGK equation in a hyperbolic scaling, which leads to the macroscopic models such as the Euler and compressible Navier-Stokes equations in the asymptotic limit. Our approaches are based on the so-called micro-macro formulation of the kinetic equation which involves a natural decomposition of the problem to the equilibrium and the non-equilibrium parts. The proposed methods are formulated for the BGK equation with constant or spatially variant Knudsen number. The new ingredients for the proposed methods to achieve high order accuracy are the following: we introduce discontinuous Galerkin (DG) discretization of arbitrary order of accuracy with nodal Lagrangian basis functions in space; we employ a high order globally stiffly accurate implicit-explicit (IMEX) Runge-Kutta (RK) scheme as time discretization. Two versions of the schemes are proposed: Scheme I is a direct formulation based on the micro-macro decomposition of the BGK equation, while Scheme II, motivated by the asymptotic analysis for the continuous problem, utilizes certain properties of the projection operator. Compared with Scheme I, Scheme II not only has better computational efficiency (the computational cost is reduced by half roughly), but also allows the establishment of a formal asymptotic analysis. Specifically, it is demonstrated that when $0<\epsilon \ll 1$, Scheme II, up to $\mathcal{O}(\epsilon^2)$, becomes a local DG discretization with an explicit RK method for the macroscopic compressible Navier-Stokes equations, a method in a similar spirit to the ones in [Bassi \& Rabey 1997, Cockburn \& Shu 1998]. Numerical results are presented for a wide range of Knudsen number to illustrate the effectiveness and high order accuracy of the methods.

In this paper, we investigate the application of the maximum principle preserving (MPP) parametrized flux limiters to the high order finite volume scheme with Runge-Kutta time discretization for solving convection dominated problems. Such flux limiter was originally proposed in [Xu, Math. Comp., 2013] and further developed in [Xiong et. al., J. Comp. Phys., 2013] for finite difference WENO schemes with Runge-Kutta time discretization for convection equations. The main idea is to limit the temporal integrated high order numerical flux toward a first order MPP monotone flux. In this paper, we generalize such flux limiter to high order finite volume methods solving convection-dominated problems, which is easy to implement and introduces little computational overhead. More importantly, for the first time in the finite volume setting, we provide a general proof that the proposed flux limiter maintains high order accuracy of the original WENO scheme for linear advection problems without any additional time step restriction. For general nonlinear convection-dominated problems, we prove that the proposed flux limiter introduces up to $O(\Delta x^3+\Delta t^3)$ modification to the high order temporal integrated flux in the original WENO scheme without extra time step constraint. We also numerically investigate the preservation of up to ninth order accuracy of the proposed flux limiter in a general setting. The advantage of the proposed method is demonstrated through various numerical experiments.

We propose a bound-preserving Runge-Kutta (RK) discontinuous Galerkin (DG) method as an efficient, effective and compact numerical approach for numerical simulation of traffic flow problems on networks, with arbitrary high order accuracy. Road networks are modeled by graphs, composed of a finite number of roads that meet at junctions. On each road, a scalar conservation law describes the dynamics, while coupling conditions are specified at junctions to define flow separation or convergence at the points where roads meet. We incorporate such coupling conditions in the RK DG framework, and apply an arbitrary high order bound preserving limiter to the RK DG method to preserve the physical bounds on the network solutions (car density). We showcase the proposed algorithm on several benchmark test cases from the literature, as well as several new challenging examples with rich solution structures. Modeling and simulation of Cauchy problems for traffic flows on networks is notorious for lack of uniqueness or (Lipschitz) continuous dependence. The discontinuous Galerkin method proposed here deals elegantly with these problems, and is perhaps the only realistic and efficient high-order method for network problems.

In this paper, we develop parametrized positivity satisfying flux limiters for the high order finite difference Runge-Kutta weighted essentially non-oscillatory (WENO) scheme solving compressible Euler equations to maintain positive density and pressure. Negative density and pressure, which often leads to simulation blow-ups or nonphysical solutions, emerges from many high resolution computations in some extreme cases. The methodology we propose in this paper is a nontrivial generalization of the parametrized maximum principle preserving flux limiters for high order finite difference schemes solving scalar hyperbolic conservation laws [22, 10, 20]. To preserve the maximum principle, the high order flux is limited towards a first order monotone flux, where the limiting procedures are designed by decoupling linear maximum principle constraints. High order schemes with such flux limiters are shown to preserve the high order accuracy via local truncation error analysis and by extensive numerical experiments with mild CFL constraints. The parametrized flux limiting approach is generalized to the Euler system to preserve the positivity of density and pressure of numerical solutions via decoupling some nonlinear constraints. Compared with existing high order positivity preserving approaches [24, 26, 25], our proposed algorithm is positivity preserving by the design; it is computationally efficient and maintains high order spatial and temporal accuracy in our extensive numerical tests. Numerical tests are performed to demonstrate the efficiency and effectiveness of the proposed new algorithm.

In this paper, we propose a Runge-Kutta (RK) central discontinuous Galerkin (CDG) gas-kinetic BGK method for the Navier-Stokes equations. The proposed method is based on the CDG method defined on two sets of overlapping meshes to avoid discontinuous solutions at cell interfaces, as well as the gas-kinetic BGK model to evaluate fluxes for both convection and diffusion terms. Redundant representation of the numerical solution in the CDG method offers great convenience in the design of gas-kinetic BGK fluxes. Specifically, the evaluation of fluxes at cell interfaces of one set of computational mesh is right inside the cells of the staggered mesh, hence the corresponding particle distribution function for flux evaluation is much simpler than that in existing gas-kinetic BGK methods. As a central scheme, the proposed CDG-BGK has doubled the memory requirement as the corresponding DG scheme; on the other hand, {for the convection part,} the CFL time step constraint of the CDG method for numerical stability is relatively large compared with that for the DG method. Numerical boundary conditions have to be treated with special care. Numerical examples for 1D and 2D viscous flow simulations are presented to validate the accuracy and robustness of the proposed RK CDG-BGK method.

In this paper, we propose the parametrized maximum principle preserving (MPP) flux limiter, originally developed in [Z. Xu, Math. Comp., (2013), in press], to the semi- Lagrangian finite difference weighted essentially non-oscillatory scheme for solving the Vlasov equation. The MPP flux limiter is proved to maintain up to fourth order accuracy for the semi-Lagrangian finite difference scheme without any time step restriction. Numerical studies on the Vlasov-Poisson system demonstrate the performance of the proposed method and its ability in preserving the positivity of the probability distribution function while maintaining the high order accuracy.

In this paper, we present error estimates of the integral deferred correction method constructed with stiffly accurate implicit Runge-Kutta methods with a nonsingular matrix $A$ in its Butcher table representation, when applied to stiff problems characterized by a small positive parameter $\varepsilon$. In our error estimates, we expand the global error in powers of $\varepsilon$ and show that the coefficients are global errors of the integral deferred correction method applied to a sequence of differential algebraic systems. A study of these errors and of the remainder of the expansion yields sharp error bounds for the stiff problem. Numerical results for the van der Pol equation are presented {to} illustrate our theoretical findings. Finally, we study the linear stability properties of these methods., Comment: 27 pages, 4 figures

Semi-Lagrangian schemes with various splitting methods, and with different reconstruction/interpolation strategies have been applied to kinetic simulations. For example, the order of spatial accuracy of the algorithms proposed in {[Qiu and Christlieb, J. Comp. Phys., 2010]} is very high (as high as ninth order). However, the temporal error is dominated by the operator splitting error, which is second order for Strang splitting. It is therefore important to overcome such low order splitting error, in order to have numerical algorithms that achieve higher orders of accuracy in both space and time. In this paper, we propose to use the integral deferred correction (IDC) method to reduce the splitting error. Specifically, the temporal order accuracy is increased by $r$ with each correction loop in the IDC framework, where $r=1,\,2$ for coupling the first order splitting and the Strang splitting, respectively. The proposed algorithm is applied to the Vlasov-Poisson system, the guiding center model, and two dimensional incompressible flow simulations in the vorticity stream-function formulation. We show numerically that the IDC procedure can automatically increase the order of accuracy in time. We also investigate numerical stability of the proposed algorithm via performing Fourier analysis to a linear model problem.