We propose a fast deterministic scheme for the space-homogeneous Boltzmann equation that exploits the low-rank structure of the velocity distribution. This paper consists of two independent contributions. The first is a lifting-projection (LP) scheme, inspired by the approach in the recent theoretical breakthroughs on the well-posedness of the Landau and Boltzmann equations. In particular, the approach lifts the nonlinear 3D Boltzmann equation to the 6D linear Kac master equation, advanced over a single time step, and projected back to its marginal in 3D. The second contribution is a low-rank tensor method for evaluating the collision operator, in which the lifted solution is represented in tensor train (TT) format and computed via a TT cross approximation algorithm with interpolation, complemented by a TT-friendly conservation correction that enforces conservation of mass, momentum, and energy. When the solution is low-rank in velocity, the method scales linearly in n when cubic interpolation is used (and quadratic in n when spectral interpolation is used), where n is the number of grid points in each velocity direction. Therefore, our methods offer significant computational savings over existing deterministic solvers in such cases. Numerical experiments on 2D and 3D benchmarks, including the BKW exact solution and anisotropic initial data, confirm the computational scaling, the expected order of accuracy and verify the effectiveness of the conservation correction.

Numerical simulation of quantum computing hardware and open quantum systems governed by the Lindblad equation is challenging due to the high dimensionality of the density matrix and the need to preserve fundamental physical properties. In our previous work, we developed an arbitrary-order, low-rank, completely positive and trace preserving (CPTP) method for the Lindblad equation with time-dependent Hamiltonians by nested Picard iteration (NPI). In this work, we develop Gregory NPI schemes, which are CPTP schemes constructed by Gregory-type quadrature on equispaced nodes. The methods, which are of order up to nine, substantially reduce the computational cost compared to our previously proposed NPI schemes with Gaussian quadrature rules, while retaining high-order accuracy and structure preservation. We analyze the stability of the resulting scheme for a physics-based test equation. Numerical experiments verify the convergence of the method and demonstrate the effectiveness of the low-rank approximation. We study the performance of a previously constructed CNOT gate for both closed and open quantum systems.

We propose a family of low-rank, completely positive and trace preserving schemes for the Lindblad equation, a common model for open quantum systems. Low-rank representation is employed at two levels: the density matrix is factorized into the product of tall-skinny matrices, and the columns of these matrices are further represented using the tensor train (TT) format, also know as matrix product states (MPS). This two-level low-rank format fits naturally into our existing Kraus is King scheme (arXiv:2409.08898v2 [math.NA]) for the Lindblad equation, whose underlying operations are arithmetic on the columns of the tall-skinny matrices. We show how these operations can be performed efficiently in the TT/MPS format, with particular emphasis on density matrix rank-truncation. We conclude with extensive numerical experiments demonstrating the convergence of this scheme and its efficiency in simulating systems with up to 10^19 degrees of freedom using only modest compute resources.

In this paper, we develop a framework for designing arbitrary high order low-rank schemes for the Lindblad equation with time-dependent Hamiltonians. Our approach is based on nested Picard iterative integrators (NPI) and results in schemes in Kraus form that are completely positive and trace preserving (CPTP). The schemes are amenable to low rank formulations, making them suitable for problems where the matrix rank of the density matrix is small.

The radiative transfer equation (RTE) is a fundamental mathematical model to describe physical phenomena involving the propagation of radiation and its interactions with the host medium, and it arises in many applications. Deterministic methods can produce accurate solutions without any statistical noise, yet often at a price of expensive computational costs originating from the intrinsic high dimensionality of the model. This is more prominent in multi-query tasks, e.g., inverse problems and optimal design, when the RTE needs to be solved repeatedly. This motivates the developments of dimensionality and model order reduction techniques for such transport models.With this work, we present the first systematic investigation of projection-based reduced order models (ROMs) following the reduced basis method (RBM) framework to simulate the parametric steady-state RTE with isotropic scattering and one energy group. The use of RBM compared to standard proper orthogonal decomposition (POD) is well motivated, especially considering that a large number of degrees of freedom is needed by full order models to solve high dimensional transport models like RTE. Four ROMs are designed, with each defining a nested family of reduced surrogate solvers of different resolution/fidelity. They are based on either a Galerkin or least-squares Petrov-Galerkin projection and utilize either an L1 or residual-based importance/error indicator. Two of the proposed ROMs are certified in the setting when the absorption cross section is positively bounded below uniformly. One technical focus and contribution lie in the proposed implementation strategies under the affine assumption of the parameter dependence of the model. These well-crafted broadly applicable strategies not only ensure the efficiency and accuracy of the offline training stage and the online prediction of reduced surrogate solvers, they also take into account the conditioning of the reduced systems as well as the stagnation-free residual evaluation for numerical robustness. Computational complexities are derived for both the offline training and online prediction stages of the proposed model order reduction strategies, and they are demonstrated numerically along with the accuracy and robustness of the reduced surrogate solvers. Numerically we observe four to six orders of magnitude speedup of our ROMs compared to full order models for some 2D2v examples.

This work proposes a new class of preconditioners for the low rank Generalized Minimal Residual Method (GMRES) for multiterm matrix equations arising from implicit timestepping of linear matrix differential equations. We are interested in computing low rank solutions to matrix equations, e.g. arising from spatial discretization of stiff partial differential equations (PDEs). The low rank GMRES method is a particular class of Krylov subspace method where the iteration is performed on the low rank factors of the solution. Such methods can exploit the low rank property of the solution to save on computational and storage cost. Of critical importance for the efficiency and applicability of the low rank GMRES method is the availability of an effective low rank preconditioner that operates directly on the low rank factors of the solution and that can limit the iteration count and the maximal Krylov rank. The preconditioner we propose here is based on the basis update and Galerkin (BUG) method, resulting from the dynamic low rank approximation. It is a nonlinear preconditioner for the low rank GMRES scheme that naturally operates on the low rank factors. Extensive numerical tests show that this new preconditioner is highly efficient in limiting iteration count and maximal Krylov rank. We show that the preconditioner performs well for general diffusion equations including highly challenging problems, e.g. high contrast, anisotropic equations. Further, it compares favorably with the state of the art exponential sum preconditioner. We also propose a hybrid BUG - exponential sum preconditioner based on alternating between the two preconditioners.

In this work, we develop implicit rank-adaptive schemes for time-dependent matrix differential equations. The dynamic low rank approximation (DLRA) is a well-known technique to capture the dynamic low rank structure based on Dirac-Frenkel time-dependent variational principle. In recent years, it has attracted a lot of attention due to its wide applicability. Our schemes are inspired by the three-step procedure used in the rank adaptive version of the unconventional robust integrator (the so called BUG integrator) for DLRA. First, a prediction (basis update) step is made computing the approximate column and row spaces at the next time level. Second, a Galerkin evolution step is invoked using a base implicit solve for the small core matrix. Finally, a truncation is made according to a prescribed error threshold. Since the DLRA is evolving the differential equation projected on to the tangent space of the low rank manifold, the error estimate of the BUG integrator contains the tangent projection (modeling) error which cannot be easily controlled by mesh refinement. This can cause convergence issue for equations with cross terms. To address this issue, we propose a simple modification, consisting of merging the row and column spaces from the explicit step truncation method together with the BUG spaces in the prediction step. In addition, we propose an adaptive strategy where the BUG spaces are only computed if the residual for the solution obtained from the prediction space by explicit step truncation method, is too large. We prove stability and estimate the local truncation error of the schemes under assumptions. We benchmark the schemes in several tests, such as anisotropic diffusion, solid body rotation and the combination of the two, to show robust convergence properties.

This paper proposes a new framework for computing low-rank solutions to nonlinear matrix equations arising from spatial discretization of nonlinear partial differential equations: low-rank Anderson acceleration (lrAA). lrAA is an adaptation of Anderson acceleration (AA), a well-known approach for solving nonlinear fixed point problems, to the low-rank format. In particular, lrAA carries out all linear and nonlinear operations in low-rank form with rank truncation using an adaptive truncation tolerance. We propose a simple scheduling strategy to update the truncation tolerance throughout the iteration according to a residual indicator. This controls the intermediate rank and iteration number effectively. To perform rank truncation for nonlinear functions, we propose a new cross approximation, which we call Cross-DEIM, with adaptive error control that is based on the discrete empirical interpolation method (DEIM). Cross-DEIM employs an iterative update between the approximate singular value decomposition (SVD) and cross approximation. It naturally incorporates a warm-start strategy for each lrAA iterate. We demonstrate the superior performance of lrAA applied to a range of linear and nonlinear problems, including those arising from finite difference discretizations of Laplace's equation, the Bratu problem, the elliptic Monge-Ampére equation and the Allen-Cahn equation.

We design high order accurate methods that exploit low rank structure in the density matrix while respecting the essential structure of the Lindblad equation. Our methods preserves complete positivity and are trace preserving.

This paper proposes two new algorithms related to the Tucker tensor format. The first method is a new cross approximation for Tucker tensors, which we call Cross^2-DEIM. Cross^2-DEIM is an iterative method that uses a fiber sampling strategy, sampling O(r) fibers in each mode, where r denotes the target rank. The fibers are selected based on the discrete empirical interpolation method (DEIM). Cross^2-DEIM resemblances the Fiber Sampling Tucker Decomposition (FSTD)2 approximation, and has favorable computational scaling compared to existing methods in the literature. We demonstrate good performance of Cross^2-DEIM in terms of iteration count and intermediate memory. First we design a fast direct Poisson solver based on Cross^2-DEIM and the fast Fourier transform. This solver can be used as a stand alone or as a preconditioner for low-rank solvers for elliptic problems. The second method is a low-rank solver for nonlinear tensor equation in Tucker format by Anderson acceleration (AA), which we call Tucker-AA. Tucker-AA is an extension of low-rank AA (lrAA) proposed in our prior work for low-rank solution to nonlinear matrix equation. We apply Cross^2-DEIM with warm-start in Tucker-AA to deal with the nonlinearity in the equation. We apply low-rank operations in AA, and by an appropriate rank truncation strategy, we are able to control the intermediate rank growth. We demonstrated the performance for Tucker-AA for approximate solutions nonlinear PDEs in 3D.

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.

This paper explores the discontinuous Galerkin (DG) methods for solving the Vlasov–Maxwell (VM) system, a fundamental model for collisionless magnetized plasma. The DG method provides an accurate numerical description with conservation and stability properties. This work studies the applicability of a post-processing technique to the DG solution in order to enhance its accuracy and resolution for the VM system. In particular, superconvergence in the negative-order norm for the probability distribution function and the electromagnetic fields is established for the DG solution. Numerical tests including Landau damping, two-stream instability, and streaming Weibel instabilities are considered showing the performance of the post-processor.

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++ software package called AdaM-DG, implementing the aSG-DG method, is available on GitHub at 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 the numerical error and the CPU cost for several benchmark tests, including linear transport equations, wave equations, and Hamilton-Jacobi (HJ) equations.

This is the third paper in a series in which we develop machine learning (ML) moment closure models for the radiative transfer equation. In our previous work (Huang et al. in J Comput Phys 453:110941, 2022), 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 (Huang et al., in: Machine learning moment closure models for the radiative transfer equation II: enforcing global hyperbolicity in gradient based closures, 2021. arXiv:2105.14410 ), 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.

In this work, we develop energy stable numerical methods to simulate electromagnetic waves propagating in optical media where the media responses include the linear Lorentz dispersion, the instantaneous nonlinear cubic Kerr response, and the nonlinear delayed Raman molecular vibrational response. Unlike the first-order PDE-ODE governing equations considered previously in Bokil et al. (J Comput Phys 350: 420–452, 2017) and Lyu et al. (J Sci Comput 89: 1–42, 2021), a model of mixed-order form is adopted here that consists of the first-order PDE part for Maxwell’s equations coupled with the second-order ODE part (i.e., the auxiliary differential equations) modeling the linear and nonlinear dispersion in the material. The main contribution is a new numerical strategy to treat the Kerr and Raman nonlinearities to achieve provable energy stability property within a second-order temporal discretization. A nodal discontinuous Galerkin (DG) method is further applied in space for efficiently handling nonlinear terms at the algebraic level, while preserving the energy stability and achieving high-order accuracy. Indeed with d_E as the number of the components of the electric field, only a d_E× d_E nonlinear algebraic system needs to be solved at each interpolation node, and more importantly, all these small nonlinear systems are completely decoupled over one time step, rendering very high parallel efficiency. We evaluate the proposed schemes by comparing them with the methods in Bokil et al. (2017) and Lyu et al. (2021) (implemented in nodal form) regarding the accuracy, computational efficiency, and energy stability, by a parallel scalability study, and also through the simulations of the soliton-like wave propagation in one dimension, as well as the spatial-soliton propagation and two-beam interactions modeled by the two-dimensional transverse electric (TE) mode of the equations.

In this paper, we propose a class of adaptive multiresolution (also called the 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 the KdV equation in [7]. For the ZK equation, which contains mixed derivative terms, we develop a new UWDG formulation. The L-2 stability is established for this new scheme on regular meshes, and the optimal error estimate with a novel local projection is obtained for a simplified ZK equation. 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.

Linear kinetic transport equations play a critical role in optical tomography, radiative transfer and neutron transport. The fundamental difficulty hampering their efficient and accurate numerical resolution lies in the high dimensionality of the physical and velocity/angular variables and the fact that the problem is multiscale in nature. Leveraging the existence of a hidden low-rank structure hinted by the diffusive limit, in this work, we design and test the angular-space reduced order model for the linear radiative transfer equation, the first such effort based on the celebrated reduced basis method (RBM). Our method is built upon a high-fidelity solver employing the discrete ordinates method in the angular space, an asymptotic preserving upwind discontinuous Galerkin method for the physical space, and an efficient synthetic accelerated source iteration for the resulting linear system. Addressing the challenge of the parameter values (or angular directions) being coupled through an integration operator, the first novel ingredient of our method is an iterative procedure where the macroscopic density is constructed from the RBM snapshots, treated explicitly and allowing a transport sweep, and then updated afterwards. A greedy algorithm can then proceed to adaptively select the representative samples in the angular space and form a surrogate solution space. The second novelty is a least-squares density reconstruction strategy, at each of the relevant physical locations, enabling the robust and accurate integration over an arbitrarily unstructured set of angular samples toward the macroscopic density. Numerical experiments indicate that our method is effective for computational cost reduction in a variety of regimes.

The Hamilton-Jacobi (HJ) equations arise in optimal control and many other applications. Oftentimes, such equations are posed in high dimensions, and this presents great numerical challenges. In this paper, we propose an adaptive sparse grid (also called adaptive multiresolution) local discontinuous Galerkin (DG) method 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 including the orthonormal Alpert & rsquo;s multiwavelets and the interpolatory multiwavelets are used to achieve multiresolution. Numerical tests in up to four dimensions are provided to validate the performance of the method. (c) 2021 Elsevier Inc. All rights reserved. Superscript/Subscript Available

This paper develops a high-order adaptive scheme for solving nonlinear Schrödinger equations. The solutions to such equations often exhibit solitary wave and local structures, which make 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.

In our recent work [Z. Peng et al., J. Comput. Phys., 415 (2020), 109485], a family of high-order asymptotic preserving (AP) methods, termed 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 subfamily 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\ll 1$ and the uniform stability property with respect to $\varepsilon$. To capture the unconditional stability, we propose a novel discrete energy and explore various stabilization mechanisms of the method and their relative contributions in different regimes. A general form of the weight function, introduced to obtain the unconditional stability for $\varepsilon\ll 1$, is also for the first time considered in such stability analysis. Based on uniform stability, a rigorous asymptotic analysis is then carried out to show the AP property.