Search

Splitting methods constitute a well-established class of numerical schemes for solving convection-diffusion-reaction problems. They have been shown to be effective in solving problems with periodic boundary conditions. However, in the case of Dirichlet boundary conditions, order reduction has been observed even with homogeneous boundary conditions. In this paper, we propose a novel splitting approach, the so-called `initial-corrected splitting method', which succeeds in overcoming order reduction. A convergence analysis is performed to demonstrate second-order convergence of this modified Strang splitting method. Furthermore, we conduct numerical experiments to illustrate the performance of the newly developed splitting approach.

In kinetic theory, numerically solving the full Boltzmann equation is extremely expensive. This is because the Boltzmann collision operator involves a high-dimensional, nonlinear integral that must be evaluated at each spatial grid point and every time step. The challenge becomes even more pronounced in the fluid (strong collisionality) regime, where the collision operator exhibits strong stiffness, causing explicit time integrators to impose severe stability restrictions. In this paper, we propose addressing this problem through a dynamical low-rank (DLR) approximation. The resulting algorithm requires evaluating the Boltzmann collision operator only $r^2$ times, where $r$, the rank of the approximation, is much smaller than the number of spatial grid points. We propose a novel DLR integrator, called the XL integrator, which reduces the number of steps compared to the available alternatives (such as the projector splitting or basis update & Galerkin (BUG) integrator). For a class of problems including the Boltzmann collision operator which enjoys a separation property between physical and velocity space, we further propose a specialized version of the XL integrator, called the sXL integrator. This version requires solving only one differential equation to update the low-rank factors. Furthermore, the proposed low-rank schemes are asymptotic-preserving, meaning they can capture the asymptotic fluid limit in the case of strong collisionality. Our numerical experiments demonstrate the efficiency and accuracy of the proposed methods across a wide range of regimes, from non-stiff (kinetic) to stiff (fluid).

Computing numerical solutions of the thermal radiative transfer equations on a finely resolved grid can be costly due to high computational and memory requirements. A numerical reduced order method that has recently been applied to a wide variety of kinetic partial differential equations is the concept of dynamical low-rank approximation (DLRA). In this paper, we consider the thermal radiative transfer equations with Su-Olson closure, leading to a linearized kinetic model. For the conducted theoretical and practical considerations we use a multiplicative splitting of the distribution function that poses additional challenges in finding an energy stable discretization and deriving a hyperbolic Courant-Friedrichs-Lewy (CFL) condition. We propose such an energy stable DLRA scheme that makes use of the augmented basis update & Galerkin integrator. This integrator allows for additional basis augmentations, enabling us to give a mathematically rigorous proof of energy stability and local mass conservation. Numerical examples confirm the derived properties and show the computational advantages of the DLRA scheme compared to a numerical solution of the full system of equations.

Boolean reaction networks are an important tool in biochemistry for studying mechanisms in the biological cell. However, the stochastic formulation of such networks requires the solution of a master equation which inherently suffers from the curse of dimensionality. In the past, the dynamical low-rank (DLR) approximation has been repeatedly used to solve high-dimensional reaction networks by separating the network into smaller partitions. However, the partitioning of these networks was so far only done by hand. In this paper, we present a heuristic, automatic partitioning scheme based on two ingredients: the Kernighan-Lin algorithm and information entropy. Our approach is computationally inexpensive and can be easily incorporated as a preprocessing step into the existing simulation workflow. We test our scheme by partitioning Boolean reaction networks on a single level and also in a hierarchical fashion with tree tensor networks. The resulting accuracy of the scheme is superior to both partitionings chosen by human experts and those found by simply minimizing the number of reaction pathways between partitions.

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.

We introduce two novel interpolatory dynamical low-rank (DLR) approximation methods for the efficient time integration of the Boltzmann-BGK equation. Both methods overcome limitations of classic DLR schemes based on orthogonal projections for nonlinear equations. In particular, we demonstrate that the proposed methods can efficiently compute solutions to the full Boltzmann-BGK equation without restricting to e.g. weakly compressible or isothermal flow. The first method we propose directly applies the recently developed interpolatory projector-splitting scheme on low-rank matrix manifolds. The second method is a variant of the rank-adaptive basis update and Galerkin scheme, where the Galerkin step is replaced by a collocation step, resulting in a new scheme we call basis update and collocate (BUC). Numerical experiments in both fluid and kinetic regimes demonstrate the performance of the proposed methods. In particular we demonstrate that the methods can be used to efficiently compute low-rank solutions in the six-dimensional (three spatial and three velocity dimensions) setting on a standard laptop.

The numerical method of dynamical low-rank approximation (DLRA) has recently been applied to various kinetic equations showing a significant reduction of the computational effort. In this paper, we apply this concept to the linear Boltzmann-Bhatnagar-Gross-Krook (Boltzmann-BGK) equation which due its high dimensionality is challenging to solve. Inspired by the special structure of the non-linear Boltzmann-BGK problem, we consider a multiplicative splitting of the distribution function. We propose a rank-adaptive DLRA scheme making use of the basis update & Galerkin integrator and combine it with an additional basis augmentation to ensure numerical stability, for which an analytical proof is given and a classical hyperbolic Courant-Friedrichs-Lewy (CFL) condition is derived. This allows for a further acceleration of computational times and a better understanding of the underlying problem in finding a suitable discretization of the system. Numerical results of a series of different test examples confirm the accuracy and efficiency of the proposed method compared to the numerical solution of the full system.

In this paper, we consider the application of exponential integrators to problems that are advection dominated, either on the entire or on a subset of the domain. In this context, we compare Leja and Krylov based methods to compute the action of exponential and related matrix functions. We set up a performance model by counting the different operations needed to implement the considered algorithms. This model assumes that the evaluation of the right-hand side is memory bound and allows us to evaluate performance in a hardware independent way. We find that exponential integrators perform comparably to explicit Runge-Kutta schemes for problems that are advection dominated in the entire domain. Moreover, they are able to outperform explicit methods in situations where small parts of the domain are diffusion dominated. We generally observe that Leja based methods outperform Krylov iterations in the problems considered. This is in particular true if computing inner products is expensive.

A constant intensity beam that propagates into a stationary plasma results in a bump-on-tail feature in velocity space. This results in an instability that transfers kinetic energy from the plasma to the electric field. We show that there are intensity profiles for the beam (found by numerical optimization) that can largely suppress this instability and drive the system into a state that, after the beam has been switched off, remains stable over long times. The modulated beam intensity requires no feedback, i.e. no knowledge of the physical system during time evolution is required, and the frequency of the modulation scales approximately inversely with system size, which is particularly favorable for large plasma systems. We also show that the results obtained are robust in the sense that perturbations, e.g. deviation from the optimized beam profiles, can be tolerated without losing the ability to suppress the instability.

In this paper we propose a semi-Lagrangian discontinuous Galerkin solver for the simulation of the scrape off layer for an electron-ion plasma. We use a time adaptive velocity space to deal with fast particles leaving the computational domain, a block structured mesh to resolve the sharp gradient in the plasma sheath, and limiters to avoid oscillations in the density function. In particular, we propose a limiter that can be computed directly from the information used in the semi-Lagrangian discontinuous Galerkin advection step. This limiter is particularly efficient on graphic processing units (GPUs) and compares favorable with limiters from the literature. We provide numerical results for a set of benchmark problems and compare different limiting strategies.

Plasma instabilities are a major concern in plasma science, for applications ranging from particle accelerators to nuclear fusion reactors. In this work, we consider the possibility of controlling such instabilities by adding an external electric field to the Vlasov--Poisson equations. Our approach to determining the external electric field is based on conducting a linear analysis of the resulting equations. We show that it is possible to select external electric fields that completely suppress the plasma instabilities present in the system when the equilibrium distribution and the perturbation are known. In fact, the proposed strategy returns the plasma to its equilibrium with a rate that is faster than exponential in time. We further perform numerical simulations of the nonlinear two-stream and bump-on-tail instabilities to verify our theory and to compare the different strategies that we propose in this work.

The stochastic description of chemical reaction networks with the kinetic chemical master equation (CME) is important for studying biological cells, but it suffers from the curse of dimensionality: The amount of data to be stored grows exponentially with the number of chemical species and thus exceeds the capacity of common computational devices for realistic problems. Therefore, time-dependent model order reduction techniques such as the dynamical low-rank approximation are desirable. In this paper we propose a dynamical low-rank algorithm for the kinetic CME using binary tree tensor networks. The dimensionality of the problem is reduced in this approach by hierarchically dividing the reaction network into partitions. Only reactions that cross partitions are subject to an approximation error. We demonstrate by two numerical examples (a 5-dimensional lambda phage model and a 20-dimensional reaction cascade) that the proposed method drastically reduces memory consumption and shows improved computational performance and better accuracy compared to a Monte Carlo method.

Dynamical low-rank approximation has become a valuable tool to perform an on-the-fly model order reduction for prohibitively large matrix differential equations. A core ingredient is the construction of integrators that are robust to the presence of small singular values and the resulting large time derivatives of the orthogonal factors in the low-rank matrix representation. Recently, the robust basis-update & Galerkin (BUG) class of integrators has been introduced. These methods require no steps that evolve the solution backward in time, often have favourable structure-preserving properties, and allow for parallel time-updates of the low-rank factors. The BUG framework is flexible enough to allow for adaptations to these and further requirements. However, the BUG methods presented so far have only first-order robust error bounds. This work proposes a second-order BUG integrator for dynamical low-rank approximation based on the midpoint rule. The integrator first performs a half-step with a first-order BUG integrator, followed by a Galerkin update with a suitably augmented basis. We prove a robust second-order error bound which in addition shows an improved dependence on the normal component of the vector field. These rigorous results are illustrated and complemented by a number of numerical experiments.

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.

Numerical simulations of kinetic problems can become prohibitively expensive due to their large memory footprint and computational costs. A method that has proven to successfully reduce these costs is the dynamical low-rank approximation (DLRA). One key question when using DLRA methods is the construction of robust time integrators that preserve the invariances and associated conservation laws of the original problem. In this work, we demonstrate that the augmented basis update & Galerkin integrator (BUG) preserves solution invariances and the associated conservation laws when using a conservative truncation step and an appropriate time and space discretization. We present numerical comparisons to existing conservative integrators and discuss advantages and disadvantages

We present an open-source CUDA-based package that consists of a compilation of exponential integrators where the action of the matrix exponential or the $\varphi_l$ functions on a vector is approximated using the method of polynomial interpolation at Leja points. Using a couple of test examples on an NVIDIA A100 GPU, we show that one can achieve significant speedups using CUDA over the corresponding CPU code. LeXInt, written in a modular format, facilitates easy integration into any existing software package, and can be used for temporal integration of any differential equation., Comment: Open-source code available at https://github.com/Pranab-JD/LeXInt; Comments & suggestions are highly welcome; if you need help using the software, please reach out

It is increasingly realized that taking stochastic effects into account is important in order to study biological cells. However, the corresponding mathematical formulation, the chemical master equation (CME), suffers from the curse of dimensionality and thus solving it directly is not feasible for most realistic problems. In this paper we propose a dynamical low-rank algorithm for the CME that reduces the dimensionality of the problem by dividing the reaction network into partitions. Only reactions that cross partitions are subject to an approximation error (everything else is computed exactly). This approach, compared to the commonly used stochastic simulation algorithm (SSA, a Monte Carlo method), has the advantage that it is completely noise-free. This is particularly important if one is interested in resolving the tails of the probability distribution. We show that in some cases (e.g. for the lambda phage) the proposed method can drastically reduce memory consumption and run time and provide better accuracy than SSA.

Anisotropic diffusion is imperative in understanding cosmic ray diffusion across the Galaxy, the heliosphere, and the interplay of cosmic rays with the Galactic magnetic field. This diffusion term contributes to the highly stiff nature of the cosmic ray transport equation. To conduct numerical simulations of time-dependent cosmic ray transport, implicit integrators (namely, Crank-Nicolson (CN)) have been traditionally favoured over the CFL-bound explicit integrators in order to be able to take large step sizes. We propose exponential methods to treat the linear anisotropc diffusion equation in the presence of advection and time-independent and time-dependent sources. These methods allow us to take even larger step sizes that can substantially speed-up the simulations whilst generating highly accurate solutions. In or subsequent work, we will use these exponential solvers in the Picard code to study anisotropic cosmic ray diffusion and we will consider additional physical processes such as continuous momentum losses and reacceleration., Comment: The 38th International Cosmic Ray Conference (ICRC2023)

Computational methods for thermal radiative transfer problems exhibit high computational costs and a prohibitive memory footprint when the spatial and directional domains are finely resolved. A strategy to reduce such computational costs is dynamical low-rank approximation (DLRA), which represents and evolves the solution on a low-rank manifold, thereby significantly decreasing computational and memory requirements. Efficient discretizations for the DLRA evolution equations need to be carefully constructed to guarantee stability while enabling mass conservation. In this work, we focus on the Su-Olson closure leading to a linearized internal energy model and derive a stable discretization through an implicit coupling of internal energy and particle density. Moreover, we propose a rank-adaptive strategy to preserve local mass conservation. Numerical results are presented which showcase the accuracy and efficiency of the proposed low-rank method compared to the solution of the full system.

We propose a dynamical low-rank algorithm for a gyrokinetic model that is used to describe strongly magnetized plasmas. The low-rank approximation is based on a decomposition into variables parallel and perpendicular to the magnetic field, as suggested by the physics of the underlying problem. We show that the resulting scheme exactly recovers the dispersion relation even with rank 1. We then perform a simulation of kinetic shear Alfv\'en waves and show that using the proposed dynamical low-rank algorithm a drastic reduction (multiple orders of magnitude) in both computational time and memory consumption can be achieved. We also compare the performance of robust first and second-order projector splitting, BUG (also called unconventional), and augmented BUG integrators as well as a FFT-based spectral and Lax--Wendroff discretization.

Fusion energy offers the potential for the generation of clean, safe, and nearly inexhaustible energy. While notable progress has been made in recent years, significant challenges persist in achieving net energy gain. Improving plasma confinement and stability stands as a crucial task in this regard and requires optimization and control of the plasma system. In this work, we deploy a PDE-constrained optimization formulation that uses a kinetic description for plasma dynamics as the constraint. This is to optimize, over all possible controllable external electric fields, the stability of the plasma dynamics under the condition that the Vlasov--Poisson (VP) equation is satisfied. For computing the functional derivative with respect to the external field in the optimization updates, the adjoint equation is derived. Furthermore, in the discrete setting, where we employ the semi-Lagrangian method as the forward solver, we also explicitly formulate the corresponding adjoint solver and the gradient as the discrete analogy to the adjoint equation and the Frechet derivative. A distinct feature we observed of this constrained optimization is the complex landscape of the objective function and the existence of numerous local minima, largely due to the hyperbolic nature of the VP system. To overcome this issue, we utilize a gradient-accelerated genetic algorithm, leveraging the advantages of the genetic algorithm's exploration feature to cover a broader search of the solution space and the fast local convergence aided by the gradient information. We show that our algorithm obtains good electric fields that are able to maintain a prescribed profile in a beam shaping problem and uses nonlinear effects to suppress plasma instability in a two-stream configuration.

In this paper we consider an approach to improve the performance of exponential Runge--Kutta integrators and Lawson schemes} in cases where the solution of a related, but usually much simpler, problem can be computed efficiently. While for implicit methods such an approach is common (e.g. by using preconditioners), for exponential integrators this has proven more challenging. Here we propose to extract a constant coefficient differential operator from the semilinear advection-diffusion-reaction equation for which, in many situations, efficient methods are known to compute the required matrix functions. Both a linear stability analysis and {\color{black} extensive} numerical experiments show that the resulting schemes can be unconditionally stable. In fact, we find that exponential integrators of Runge--Kutta type and Lawson schemes can have better stability properties than similarly constructed implicit-explicit schemes. We also derive two new Lawson type integrators that further improve on these stability properties. The overall effectiveness of the approach is highlighted by a number of performance comparisons on examples in two and three space dimensions.

Radiation transport problems are posed in a high-dimensional phase space, limiting the use of finely resolved numerical simulations. An emerging tool to efficiently reduce computational costs and memory footprint in such settings is dynamical low-rank approximation (DLRA). Despite its efficiency, numerical methods for DLRA need to be carefully constructed to guarantee stability while preserving crucial properties of the original problem. Important physical effects that one likes to preserve with DLRA include capturing the diffusion limit in the high-scattering regimes as well as dissipating energy. In this work we propose and analyze a dynamical low-rank method based on the "unconventional" basis-update & Galerkin step integrator. We show that this method is asymptotic-preserving, i.e., it captures the diffusion limit, and energy stable under a CFL condition. The derived CFL condition captures the transition from the hyperbolic to the parabolic regime when approaching the diffusion limit.

In this paper, we demonstrate the efficiency of using semi-Lagrangian discontinuous Galerkin methods to solve the drift-kinetic equation using graphic processing units (GPUs). In this setting we propose a second order splitting scheme and a 2d semi-Lagrangian scheme in the poloidal plane. The resulting method is able to conserve mass up to machine precision, allows us to take large time steps due to the absence of a CFL condition and provides local data dependency which is essential to obtain good performance on state-of-the art high-performance computing systems. We report simulations of a drift-kinetic ion temperature gradient (ITG) instability and show that our implementation achieves a performance of up to 600 GB/s on an A100 GPU.

The anisotropic diffusion equation is imperative in understanding cosmic ray diffusion across the Galaxy, the heliosphere, and its interplay with the ambient magnetic field. This diffusion term contributes to the highly stiff nature of the CR transport equation. In order to conduct numerical simulations of time-dependent cosmic ray transport, implicit integrators have been traditionally favoured over the CFL-bound explicit integrators in order to be able to take large step sizes. We propose exponential methods that directly compute the exponential of the matrix to solve the linear anisotropic diffusion equation. These methods allow us to take even larger step sizes; in certain cases, we are able to choose a step size as large as the simulation time, i.e., only one time step. This can substantially speed-up the simulations whilst generating highly accurate solutions (l2 error $\leq 10^{-10}$). Additionally, we test an approach based on extracting a constant diffusion coefficient from the anisotropic diffusion equation, where the constant coefficient term is solved implicitly or exponentially and the remainder is treated using some explicit method. We find that this approach, for homogeneous linear problems, is unable to improve on the exponential-based methods that directly evaluate the matrix exponential., Comment: in submission

Krylov-based algorithms have long been preferred to compute the matrix exponential and exponential-like functions appearing in exponential integrators. Of late, direct polynomial interpolation of the action of these exponential-like functions have been shown to be competitive with the Krylov methods. We analyse the performance of the state-of-the-art Krylov algorithm, KIOPS, and the method of polynomial interpolation at Leja points for a number of exponential integrators for various test problems and with varying amounts of stiffness. Additionally, we investigate the performance of an iterative scheme that combines both the KIOPS and Leja approach, named LeKry, that shows substantial improvements over both the Leja- and Krylov-based methods for certain exponential integrators. Whilst we do manage to single out a favoured iterative scheme for each of the exponential integrators that we consider in this study, we do not find any conclusive evidence for preferring either KIOPS or Leja for different classes of exponential integrators. We are unable to identify a superior exponential integrator, one that performs better than all others, for most, if not all of the problems under consideration. We, however, do find that the performance significantly depends on the interplay between the iterative scheme and the specific exponential integrator under consideration., Comment: in submission; attached supplemental material

We present a publicly available software for exponential integrators that computes the $\varphi_l(z)$ functions using polynomial interpolation. The interpolation method at Leja points have recently been shown to be competitive with the traditionally-used Krylov subspace method. The developed framework facilitates easy adaptation into any Python software package for time integration., Comment: Publicly available software available at https://github.com/Pranab-JD/LeXInt, in submission

Dynamical low-rank approximation, as has been demonstrated recently, can be extremely efficient in solving kinetic equations. However, a major deficiency is that they do not preserve the structure of the underlying physical problem. For example, the classic dynamical low-rank methods violate mass, momentum, and energy conservation. In [L. Einkemmer, I. Joseph, J. Comput. Phys. 443:110495, 2021] a conservative dynamical low-rank approach has been proposed. However, directly integrating the resulting equations of motion, similar to the classic dynamical low-rank approach, results in an ill-posed scheme. In this work we propose a robust, i.e. well-posed, integrator for the conservative dynamical low-rank approach that conserves mass and momentum (up to machine precision) and significantly improves energy conservation. We also report improved qualitative results for some problems and show how the approach can be combined with a rank adaptive scheme.

Running kinetic plasma physics simulations using grid-based solvers is very demanding both in terms of memory as well as computational cost. This is primarily due to the up to six-dimensional phase space and the associated unfavorable scaling of the computational cost as a function of grid spacing (often termed the curse of dimensionality). In this paper, we present 4d, 5d, and 6d simulations of the Vlasov--Poisson equation with a split-step semi-Lagrangian discontinuous Galerkin scheme on graphic processing units (GPUs). The local communication pattern of this method allows an efficient implementation on large-scale GPU-based systems and emphasizes the importance of considering algorithmic and high-performance computing aspects in unison. We demonstrate a single node performance above 2 TB/s effective memory bandwidth (on a node with 4 A100 GPUs) and show excellent scaling (parallel efficiency between 30% and 67%) for up to 1536 A100 GPUs on JUWELS Booster., Comment: Submitted to The International Journal of High Performance Computing Applications

Running kinetic simulations using grid-based methods is extremely expensive due to the up to six-dimensional phase space. Recently, it has been shown that dynamical low-rank algorithms can drastically reduce the required computational effort, while still accurately resolving important physical features such as filamentation and Landau damping. In this paper, we propose a new second order projector-splitting dynamical low-rank algorithm for the full six-dimensional Vlasov--Poisson equations. An exponential integrator based Fourier spectral method is employed to obtain a numerical scheme that is CFL condition free but still fully explicit. The resulting method is implemented with the aid of Ensign, a software framework which facilitates the efficient implementation of dynamical low-rank algorithms on modern multi-core CPU as well as GPU based systems. Its usage and features are briefly described in the paper as well. The presented numerical results demonstrate that 6D simulations can be run on a single workstation and highlight the significant speedup that can be obtained using GPUs.

We propose a novel algorithm for the temporal integration of the resistive magnetohydrodynamics (MHD) equations. The approach is based on exponential Rosenbrock schemes in combination with Leja interpolation. It naturally preserves Gauss's law for magnetism and is unencumbered by the stability constraints observed for explicit methods. Remarkable progress has been achieved in designing exponential integrators and computing the required matrix functions efficiently. However, employing them in MHD simulations of realistic physical scenarios requires a matrix-free implementation. We show how an efficient algorithm based on Leja interpolation that only uses the right-hand side of the differential equation (i.e. matrix free), can be constructed. We further demonstrate that it outperforms Krylov-based exponential integrators as well as explicit and implicit methods using test models of magnetic reconnection and the Kelvin--Helmholtz instability. Furthermore, an adaptive step-size strategy that gives excellent and predictable performance, particularly in the lenient- to intermediate-tolerance regime that is often of importance in practical applications, is employed., Comment: Accepted for publication in ApJS; 20 pages, 13 figures

The dynamical low-rank approximation (DLRA) is used to treat high-dimensional problems that arise in such diverse fields as kinetic transport and uncertainty quantification. Even though it is well known that certain spatial and temporal discretizations when combined with the DLRA approach can result in numerical instability, this phenomenon is poorly understood. In this paper we perform a $L^2$ stability analysis for the corresponding nonlinear equations of motion. This reveals the source of the instability for the projector splitting integrator when first discretizing the equations and then applying the DLRA. Based on this we propose a projector splitting integrator, based on applying DLRA to the continuous system before performing the discretization, that recovers the classic CFL condition. We also show that the unconventional integrator has more favorable stability properties and explain why the projector splitting integrator performs better when approximating higher moments, while the unconventional integrator is generally superior for first order moments. Furthermore, an efficient and stable dynamical low-rank update for the scattering term in kinetic transport is proposed. Numerical experiments for kinetic transport and uncertainty quantification, which confirm the results of the stability analysis, are presented.

Quantifying uncertainties in hyperbolic equations is a source of several challenges. First, the solution forms shocks leading to oscillatory behaviour in the numerical approximation of the solution. Second, the number of unknowns required for an effective discretization of the solution grows exponentially with the dimension of the uncertainties, yielding high computational costs and large memory requirements. An efficient representation of the solution via adequate basis functions permits to tackle these difficulties. The generalized polynomial chaos (gPC) polynomials allow such an efficient representation when the distribution of the uncertainties is known. These distributions are usually only available for input uncertainties such as initial conditions, therefore the efficiency of this ansatz can get lost during runtime. In this paper, we make use of the dynamical low-rank approximation (DLRA) to obtain a memory-wise efficient solution approximation on a lower dimensional manifold. We investigate the use of the matrix projector-splitting integrator and the unconventional integrator for dynamical low-rank approximation, deriving separate time evolution equations for the spatial and uncertain basis functions, respectively. This guarantees an efficient approximation of the solution even if the underlying probability distributions change over time. Furthermore, filters to mitigate the appearance of spurious oscillations are implemented, and a strategy to enforce boundary conditions is introduced. The proposed methodology is analyzed for Burgers' equation equipped with uncertain initial values represented by a two-dimensional random vector. The numerical results show a reduction of the memory requirements, and that the important characteristics of the original system are well captured.

Recently a splitting approach has been presented for the simulation of sonic-boom propagation. Splitting methods allow one to divide complicated partial differential equations into simpler parts that are solved by specifically tailored numerical schemes. The present work proposes a second order exponential integrator for the numerical solution of sonic-boom propagation modelled through a dispersive equation with Burgers' nonlinearity. The linear terms are efficiently solved in frequency space through FFT, while the nonlinear terms are efficiently solved by a WENO scheme. The numerical method is designed to be highly parallelisable and therefore takes full advantage of modern computer hardware. The new approach also improves the accuracy compared to the splitting method and it reduces oscillations. The enclosed numerical results illustrate that parallelisation on a CPU results in a speedup of 22 times faster than the straightforward sequential version. The GPU implementation further accelerates the runtime by a factor 3, which improves to 5 when single precision is used instead of double precision.

In this paper, we propose a $\mu$-mode integrator for computing the solution of stiff evolution equations. The integrator is based on a $d$-dimensional splitting approach and uses exact (usually precomputed) one-dimensional matrix exponentials. We show that the action of the exponentials, i.e. the corresponding batched matrix-vector products, can be implemented efficiently on modern computer systems. We further explain how $\mu$-mode products can be used to compute spectral transforms efficiently even if no fast transform is available. We illustrate the performance of the new integrator by solving, among the others, three-dimensional linear and nonlinear Schr\"odinger equations, and we show that the $\mu$-mode integrator can significantly outperform numerical methods well established in the field. We also discuss how to efficiently implement this integrator on both multi-core CPUs and GPUs. Finally, the numerical experiments show that using GPUs results in performance improvements between a factor of $10$ and $20$, depending on the problem.

Traditional step size controllers make the tacit assumption that the cost of a time step is independent of the step size. This is reasonable for explicit and implicit integrators that use direct solvers. In the context of exponential integrators, however, an iterative approach, such as Krylov methods or polynomial interpolation, to compute the action of the required matrix functions is usually employed. In this case, the assumption of constant cost is not valid. This is, in particular, a problem for higher-order exponential integrators, which are able to take relatively large time steps based on accuracy considerations. In this paper, we consider an adaptive step size controller for exponential Rosenbrock methods that determines the step size based on the premise of minimizing computational cost. The largest allowed step size, given by accuracy considerations, merely acts as a constraint. We test this approach on a range of nonlinear partial differential equations. Our results show significant improvements (up to a factor of 4 reduction in the computational cost) over the traditional step size controller for a wide range of tolerances., Comment: 23 pages, 14 figures, matches with the published version

The primary challenge in solving kinetic equations, such as the Vlasov equation, is the high-dimensional phase space. In this context, dynamical low-rank approximations have emerged as a promising way to reduce the high computational cost imposed by such problems. However, a major disadvantage of this approach is that the physical structure of the underlying problem is not preserved. In this paper, we propose a dynamical low-rank algorithm that conserves mass, momentum, and energy as well as the corresponding continuity equations. We also show how this approach can be combined with a conservative time and space discretization.

It has recently been demonstrated that dynamical low-rank algorithms can provide robust and efficient approximation to a range of kinetic equations. This is true especially if the solution is close to some asymptotic limit where it is known that the solution is low-rank. A particularly interesting case is the fluid dynamic limit that is commonly obtained in the limit of small Knudsen number. However, in this case the Maxwellian which describes the corresponding equilibrium distribution is not necessarily low-rank; because of this, the methods known in the literature are only applicable to the weakly compressible case. In this paper, we propose an efficient dynamical low-rank integrator that can capture the fluid limit -- the Navier-Stokes equations -- of the Boltzmann-BGK model even in the compressible regime. This is accomplished by writing the solution as $f=Mg$, where $M$ is the Maxwellian and the low-rank approximation is only applied to $g$. To efficiently implement this decomposition within a low-rank framework requires, in the isothermal case, that certain coefficients are evaluated using convolutions, for which fast algorithms are known. Using the proposed decomposition also has the advantage that the rank required to obtain accurate results is significantly reduced compared to the previous state of the art. We demonstrate this by performing a number of numerical experiments and also show that our method is able to capture sharp gradients/shock waves.

Rational exponential integrators (REXI) are a class of numerical methods that are well suited for the time integration of linear partial differential equations with imaginary eigenvalues. Since these methods can be parallelized in time (in addition to the spatial parallelization that is commonly performed) they are well suited to exploit modern high performance computing systems. In this paper, we propose a novel REXI scheme that drastically improves accuracy and efficiency. The chosen approach will also allow us to easily determine how many terms are required in the approximation in order to obtain accurate results. We provide comparative numerical simulations for a shallow water equation that highlight the efficiency of our approach and demonstrate that REXI schemes can be efficiently implemented on graphic processing units., Comment: 25 pages, 5 figures, submitted to Journal of Computational Physics

The present work proposes a second-order time splitting scheme for a linear dispersive equation with a variable advection coefficient subject to transparent boundary conditions. For its spatial discretization, a dual Petrov--Galerkin method is considered which gives spectral accuracy. The main difficulty in constructing a second-order splitting scheme in such a situation lies in the compatibility condition at the boundaries of the sub-problems. In particular, the presence of an inflow boundary condition in the advection part results in order reduction. To overcome this issue a modified Strang splitting scheme is introduced that retains second-order accuracy. For this numerical scheme a stability analysis is conducted. In addition, numerical results are shown to support the theoretical derivations.

We introduce a dynamical low-rank method to reduce the computational complexity for solving the multi-scale multi-dimensional linear transport equation. The method is based on a macro-micro decomposition of the equation. The proposed numerical method uses the low rank approximation only for the micro part of the solution. The time and spatial discretizations are done properly so that the overall scheme is second order accurate and asymptotic-preserving (AP); that is, in the diffusive regime, the scheme becomes a macroscopic solver for the limiting diffusion equation and is automatically low rank. We demonstrate the accuracy and efficiency of the proposed low rank method by a number of two-dimensional examples.

The UPC programming language offers parallelism via logically partitioned shared memory, which typically spans physically disjoint memory sub-systems. One convenient feature of UPC is its ability to automatically execute between-thread data movement, such that the entire content of a shared data array appears to be freely accessible by all the threads. The programmer friendliness, however, can come at the cost of substantial performance penalties. This is especially true when indirectly indexing the elements of a shared array, for which the induced between-thread data communication can be irregular and have a fine-grained pattern. In this paper we study performance enhancement strategies specifically targeting such fine-grained irregular communication in UPC. Starting from explicit thread privatization, continuing with block-wise communication, and arriving at message condensing and consolidation, we obtained considerable performance improvement of UPC programs that originally require fine-grained irregular communication. Besides the performance enhancement strategies, the main contribution of the present paper is to propose performance models for the different scenarios, in form of quantifiable formulas that hinge on the actual volumes of various data movements plus a small number of easily obtainable hardware characteristic parameters. These performance models help to verify the enhancements obtained, while also providing insightful predictions of similar parallel implementations, not limited to UPC, that also involve between-thread or between-process irregular communication. As a further validation, we also apply our performance modeling methodology and hardware characteristic parameters to an existing UPC code for solving a 2D heat equation on a uniform mesh.

The efficient numerical solution of many kinetic models in plasma physics is impeded by the stiffness of these systems. Exponential integrators are attractive in this context as they remove the CFL condition induced by the linear part of the system, which in practice is often the most stringent stability constraint. In the literature, these schemes have been found to perform well, e.g., for drift-kinetic problems. Despite their overall efficiency and their many favorable properties, most of the commonly used exponential integrators behave rather erratically in terms of the allowed time step size in some situations. This severely limits their utility and robustness.Our goal in this paper is to explain the observed behavior and suggest exponential methods that do not suffer from the stated deficiencies. To accomplish this we study the stability of exponential integrators for a linearized problem. This analysis shows that classic exponential integrators exhibit severe deficiencies in that regard. Based on the analysis conducted we propose to use Lawson methods, which can be shown not to suffer from the same stability issues. We confirm these results and demonstrate the efficiency of Lawson methods by performing numerical simulations for both the Vlasov-Poisson system and a drift-kinetic model of a ion temperature gradient instability.

In this paper, our goal is to efficiently solve the Vlasov equation on GPUs. A semi-Lagrangian discontinuous Galerkin scheme is used for the discretization. Such kinetic computations are extremely expensive due to the high-dimensional phase space. The SLDG code, which is publicly available under the MIT license abstracts the number of dimensions and uses a shared codebase for both GPU and CPU based simulations. We investigate the performance of the implementation on a range of both Tesla (V100, Titan V, K80) and consumer (GTX 1080 Ti) GPUs. Our implementation is typically able to achieve a performance of approximately 470 GB/s on a single GPU and 1600 GB/s on four V100 GPUs connected via NVLink. This results in a speedup of about a factor of ten (comparing a single GPU with a dual socket Intel Xeon Gold node) and approximately a factor of 35 (comparing a single node with and without GPUs). In addition, we investigate the effect of single precision computation on the performance of the SLDG code and demonstrate that a template based dimension independent implementation can achieve good performance regardless of the dimensionality of the problem.

Dynamical low-rank algorithms are a class of numerical methods that compute low-rank approximations of dynamical systems. This is accomplished by projecting the dynamics onto a low-dimensional manifold and writing the solution directly in terms of the low-rank factors. The approach has been successfully applied to many types of differential equations. Recently, efficient dynamical low-rank algorithms have been applied to treat kinetic equations, including the Vlasov--Poisson and the Boltzmann equation, where it was demonstrated that the methods are able to capture the low-rank structure of the solution and significantly reduce numerical cost, while often maintaining high accuracy. However, no numerical analysis is currently available. In this paper, we investigate the error analysis for a dynamical low-rank algorithm applied to the multi-scale linear Boltzmann equation (a classical model in kinetic theory) to showcase the validity of the application of dynamical low-rank algorithms to kinetic theory. The equation, in its parabolic regime, is known to be rank one theoretically, and we will prove that the scheme can dynamically and automatically capture this low-rank structure. This work thus serves as the first mathematical error analysis for a dynamical low-rank approximation applied to a kinetic problem.

The goal of the present work is to solve a linear dispersive equation with variable coefficient advection on an unbounded domain. In this setting, transparent boundary conditions are vital to allow waves to leave (or even re-enter) the, necessarily finite, computational domain. To obtain an efficient numerical scheme we discretize space using a spectral method. This allows us to drastically reduce the number of grid points required for a given accuracy. Applying a fully implicit time integrator, however, would require us to invert full matrices. This is addressed by performing an operator splitting scheme and only treating the third order differential operator, stemming from the dispersive part, implicitly; this approach can also be interpreted as an implicit-explicit scheme. However, the fact that the transparent boundary conditions are non-homogeneous and depend implicitly on the numerical solution presents a significant obstacle for the splitting/pseudo-spectral approach investigated here. We show how to overcome these difficulties and demonstrate the proposed numerical scheme by performing a number of numerical simulations.