Andrew Christlieb

David Lawlor, Yang Wang, Andrew Christlieb

We present a new deterministic algorithm for the sparse Fourier transform problem, in which we seek to identify k << N significant Fourier coefficients from a signal of bandwidth N. Previous deterministic algorithms exhibit quadratic runtime scaling, while our algorithm scales linearly with k in the average case. Underlying our algorithm are a few simple observations relating the Fourier coefficients of time-shifted samples to unshifted samples of the input function. This allows us to detect when aliasing between two or more frequencies has occurred, as well as to determine the value of unaliased frequencies. We show that empirically our algorithm is orders of magnitude faster than competing algorithms.

Andrew Christlieb, Wei Guo, Yan Jiang

In this paper, we present a novel numerical scheme for solving a class of nonlinear degenerate parabolic equations with non-smooth solutions. The proposed method relies on a special kernel based formulation of the solutions found in our early work on the method of lines transpose and successive convolution. In such a framework, a high order weighted essentially non-oscillatory (WENO) methodology and a nonlinear filter are further employed to avoid spurious oscillations. High order accuracy in time is realized by using the high order explicit strong-stability-preserving (SSP) Runge-Kutta method. Moreover, theoretical investigations of the kernel based formulation combined with an explicit SSP method indicates that the combined scheme is unconditionally stable and up to third order accuracy. Evaluation of the kernel based approach is done with a fast $\mathcal{O}(N)$ summation algorithm. The new method allows for much larger time step evolution compared with other explicit schemes with the same order accuracy, leading to remarkable computational savings.

Andrew J. Christlieb, Yuan Liu, Qi Tang et al.

In this paper, we utilize the maximum-principle-preserving flux limiting technique, originally designed for high order weighted essentially non-oscillatory (WENO) methods for scalar hyperbolic conservation laws, to develop a class of high order positivity-preserving finite difference WENO methods for the ideal magnetohydrodynamic (MHD) equations. Our schemes, under the constrained transport (CT) framework, can achieve high order accuracy, a discrete divergence-free condition and positivity of the numerical solution simultaneously. Numerical examples in 1D, 2D and 3D are provided to demonstrate the performance of the proposed method.

Andrew J. Christlieb, Xiao Feng, Yan Jiang et al.

A high-order finite difference numerical scheme is developed for the ideal magnetohydrodynamic equations based on an alternative flux formulation of the weighted essentially non-oscillatory (WENO) scheme. It computes a high-order numerical flux by a Taylor expansion in space, with the lowest-order term solved from a Riemann solver and the higher-order terms constructed from physical fluxes by limited central differences. The scheme coupled with several Riemann solvers, including a Lax-Friedrichs solver and HLL-type solvers, is developed on general curvilinear meshes in two dimensions and verified on a number of benchmark problems. In particular, a HLLD solver on Cartesian meshes is extended to curvilinear meshes with proper modifications. A numerical boundary condition for the perfect electrical conductor (PEC) boundary is derived for general geometry and verified through a bow shock flow. Numerical results also confirm the advantages of using low dissipative Riemann solvers in the current framework.

Matthew F. Causley, Hana Cho, Andrew J. Christlieb et al.

We present a new solver for nonlinear parabolic problems that is L-stable and achieves high order accuracy in space and time. The solver is built by first constructing a single-dimensional heat equation solver that uses fast O(N) convolution. This fundamental solver has arbitrary order of accuracy in space, and is based on the use of the Green's function to invert a modified Helmholtz equation. Higher orders of accuracy in time are then constructed through a novel technique known as successive convolution (or resolvent expansions). These resolvent expansions facilitate our proofs of stability and convergence, and permit us to construct schemes that have provable stiff decay. The multi-dimensional solver is built by repeated application of dimensionally split independent fundamental solvers. Finally, we solve nonlinear parabolic problems by using the integrating factor method, where we apply the basic scheme to invert linear terms (that look like a heat equation), and make use of Hermite-Birkhoff interpolants to integrate the remaining nonlinear terms. Our solver is applied to several linear and nonlinear equations including heat, Allen-Cahn, and the Fitzhugh-Nagumo system of equations in one and two dimensions.

David C. Seal, Qi Tang, Zhengfu Xu et al.

In this work we construct a high-order, single-stage, single-step positivity-preserving method for the compressible Euler equations. Space is discretized with the finite difference weighted essentially non-oscillatory (WENO) method. Time is discretized through a Lax-Wendroff procedure that is constructed from the Picard integral formulation (PIF) of the partial differential equation. The method can be viewed as a modified flux approach, where a linear combination of a low- and high-order flux defines the numerical flux used for a single-step update. The coefficients of the linear combination are constructed by solving a simple optimization problem at each time step. The high-order flux itself is constructed through the use of Taylor series and the Cauchy-Kowalewski procedure that incorporates higher-order terms. Numerical results in one- and two-dimensions are presented.

Matthew Causley, Hana Cho, Andrew Christlieb

In this work, we develop an $\mathcal{O}(N)$ implicit real space method in 1D and 2D for the Cahn Hilliard (CH) and vector Cahn Hilliard (VCH) equations, based on the Method Of Lines Transpose (MOL$^\text{T}$) formulation. This formulation results in a semi-discrete time stepping algorithm, which we prove is gradient stable in the $H^{-1}$ norm. The spatial discretization follows from dimensional splitting, and an $\mathcal{O}(N)$ matrix-free solver, which applies fast convolution to the modified Helmholtz equation. We propose a novel factorization technique, in which fourth order spatial derivatives are incorporated into the solver. The splitting error is included in the nonlinear fixed point iteration, resulting in a high order, logically Cartesian (line-by-line) update. Our method is fast, but not restricted to periodic boundaries like the fast Fourier transform (FFT). The basic solver is implemented using the Backward Euler formulation, and we extend this to both backward difference (BDF) stencils, implicit Runge Kutta (SDIRK) and spectral deferred correction (SDC) frameworks to achieve high orders of temporal accuracy. We demonstrate with numerical results that the CH, and VCH equations maintain gradient stability in one and two spatial dimensions. We also explore time-adaptivity, so that meta-stable states and ripening events can be simulated both quickly and efficiently.

Yingda Cheng, Andrew J. Christlieb, Wei Guo et al.

In plasma simulations, where the speed of light divided by a characteristic length is at a much higher frequency than other relevant parameters in the underlying system, such as the plasma frequency, implicit methods begin to play an important role in generating efficient solutions in these multi-scale problems. Under conditions of scale separation, one can rescale Maxwell's equations in such a way as to give a magneto static limit known as the Darwin approximation of electromagnetics. In this work, we present a new approach to solve Maxwell's equations based on a Method of Lines Transpose (MOL$^T$) formulation, combined with a fast summation method with computational complexity $O(N\log{N})$, where $N$ is the number of grid points (particles). Under appropriate scaling, we show that the proposed schemes result in asymptotic preserving methods that can recover the Darwin limit of electrodynamics.

Benjamin Ong, Andrew Melfi, Andrew Christlieb

In this paper, we further develop a family of parallel time integrators known as Revisionist Integral Deferred Correction methods (RIDC) to allow for the semi-implicit solution of time dependent PDEs. Additionally, we show that our semi-implicit RIDC algorithm can harness the computational potential of multiple general purpose graphical processing units (GPUs) in a single node by utilizing existing CUBLAS libraries for matrix linear algebra routines in our implementation. In the numerical experiments, we show that our implementation computes a fourth order solution using four GPUs and four CPUs in approximately the same wall clock time as a first order solution computed using a single GPU and a single CPU.

Matthew Causley, Andrew Christlieb, Eric Wolf

Building upon recent results obtained in [7,8,9], we describe an efficient second order, A-stable scheme for solving the wave equation, based on the method of lines transpose (MOL$^T$), and the resulting semi-discrete (i.e. continuous in space) boundary value problem. In [7], A-stable schemes of high order were derived, and in [9] a high order, fast $\mathcal{O}(N)$ spatial solver was derived, which is matrix-free and is based on dimensional-splitting. In this work, are interested in building a wave solver, and our main concern is the development of boundary conditions. We demonstrate all desired boundary conditions for a wave solver, including outflow boundary conditions, in 1D and 2D. The scheme works in a logically Cartesian fashion, and the boundary points are embedded into the regular mesh, without incurring stability restrictions, so that boundary conditions are imposed without any reduction in the order of accuracy. We demonstrate how the embedded boundary approach works in the cases of Dirichlet and Neumann boundary conditions. Further, we develop outflow and periodic boundary conditions for the MOL$^T$ formulation. Our solver is designed to couple with particle codes, and so special attention is also paid to the implementation of point sources, and soft sources which can be used to launch waves into waveguides.

David Bortz, Andrew Christlieb

We propose a novel approach which employs random sampling to generate an accurate non-uniform mesh for numerically solving Partial Differential Equation Boundary Value Problems (PDE-BVP's). From a uniform probability distribution U over a 1D domain, we sample M discretizations of size N where M>>N. The statistical moments of the solutions to a given BVP on each of the M ultra-sparse meshes provide insight into identifying highly accurate non-uniform meshes. Essentially, we use the pointwise mean and variance of the coarse-grid solutions to construct a mapping Q(x) from uniformly to non-uniformly spaced mesh-points. The error convergence properties of the approximate solution to the PDE-BVP on the non-uniform mesh are superior to a uniform mesh for a certain class of BVP's. In particular, the method works well for BVP's with locally non-smooth solutions. We present a framework for studying the sampled sparse-mesh solutions and provide numerical evidence for the utility of this approach as applied to a set of example BVP's. We conclude with a discussion of how the near-perfect paralellizability of our approach suggests that these strategies have the potential for highly efficient utilization of massively parallel multi-core technologies such as General Purpose Graphics Processing Units (GPGPU's). We believe that the proposed algorithm is beyond embarrassingly parallel; implementing it on anything but a massively multi-core architecture would be scandalous.

Mandela B. Quashie, J. W. Burby, Andrew J. Christlieb et al.

We revisit the Scovel-Weinstein framework (Scovel & Weinstein, CPAM 1994) for reducing the Vlasov-Poisson system while preserving its Hamiltonian structure. Standard particle-in-cell (PIC) algorithms approximate the distribution function by macro-particles with position and velocity. In contrast, Scovel-Weinstein decorated particles involve additional shape degrees of freedom, while maintaining a finite-dimensional reduction with Hamiltonian structure inherited from the continuum model. Although the original work established this structure three decades ago, its computational potential has remained largely unexplored. We present a practical implementation of the Scovel-Weinstein model and compare it with a standard PIC algorithm. Numerical experiments demonstrate that macro-particles in standard PIC can be replaced by far fewer decorated particles while retaining comparable accuracy. This decorated particle approach offers a new structure-preserving paradigm for kinetic plasma simulation.

David C. Seal, Yaman Güçlü, Andrew J. Christlieb

High-order temporal discretizations for hyperbolic conservation laws have historically been formulated as either a method of lines (MOL) or a Lax-Wendroff method. In the MOL viewpoint, the partial differential equation is treated as a large system of ordinary differential equations (ODEs), where an ODE tailored time-integrator is applied. In contrast, Lax-Wendroff discretizations immediately convert Taylor series in time to discrete spatial derivatives. In this work, we propose the Picard integral formulation (PIF), which is based on the method of modified fluxes, and is used to derive new Taylor and Runge-Kutta (RK) methods. In particular, we construct a new class of conservative finite difference methods by applying WENO reconstructions to the so-called "time-averaged" fluxes. Our schemes are automatically conservative under any modification of the fluxes, which is attributed to the fact that classical WENO reconstructions conserve mass when coupled with forward Euler time steps. The proposed Lax-Wendroff discretization is constructed by taking Taylor series of the flux function as opposed to Taylor series of the conserved variables. The RK discretization differs from classical MOL formulations because we apply WENO reconstructions to time-averaged fluxes rather than taking linear combinations of spatial derivatives of the flux. In both cases, we only need one projection onto the characteristic variables per time step. The PIF is generic, and lends itself to a multitude of options for further investigation. At present, we present two canonical examples: one based on Taylor, and the other based on the classical RK method. Stability analyses are presented for each method. The proposed schemes are applied to hyperbolic conservation laws in one- and two-dimensions and the results are in good agreement with current state of the art methods.

Andrew J. Christlieb, Sining Gong, Jing-Mei Qiu et al.

The Wigner-Poisson system is a deterministic phase-space model for quantum kinetic electron dynamics, but high-dimensional simulations are limited by the full 3D3V phase space and the nonlocal Wigner potential. We develop a structure-preserving, sampling-based adaptive-rank solver in hierarchical Tucker format for finite-$H$ regimes in which Wigner-Poisson solutions exhibit exploitable low-rank structure. The central difficulty is that adaptive compression can destroy the Fourier-Hermitian tensor symmetry required for a real inverse velocity transform and can break discrete global conservation laws. We address these issues with a Fourier-Hermitian-symmetry-aware sampling and mapping procedure and a global moment correction enforcing mass, momentum, and self-consistent total energy. Numerical tests for two-stream instability and strong Landau damping in 2D2V and 3D3V show roundoff-level conservation, preservation of the real-valued inverse transform, and approximately linear scaling with respect to the number of grid points per coordinate over the tested rank range. The results demonstrate that long-time 3D3V Wigner-Poisson simulations can be performed without assembling the full phase-space tensor.

Andrew J. Christlieb, Mingchang Ding, Juntao Huang et al.

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

Andrew J. Christlieb, Luis Chacon, Sining Gong

We develop an unstaggered, potential-based particle-in-cell method for the nonrelativistic Vlasov-Maxwell system in the Lorenz gauge. The field update is written as a Crank-Nicolson discretization of first-order wave systems for the scalar potential, the vector potential, and their time derivatives. The charge density is not deposited directly; instead, it is advanced from the discrete continuity equation using the current deposited from the particles. This opens up algorithmic flexibility with a range of innovation, including unstaggered mesh layouts that preserve the Lorenz gauge and Gauss's law at the discrete level. In the potential formulation, this source ordering also permits preservation of the Lorenz gauge and Gauss's law at the discrete level. To extend the paradigm to an energy-conserving formulation, we introduce a consistent orbit-averaged scatter, gather, and particle push. For energy consistency, the update of the canonical momentum is modified by replacing the pointwise midpoint derivative of the vector potential with an orbit-averaged discrete gradient of the mesh-interpolated vector potential consistent with the orbit-average maps. This construction satisfies an exact finite-difference chain rule along each particle orbit. As a result, the particle work equals the mesh work appearing in the Crank-Nicolson field-energy balance, yielding exact total-energy conservation up to nonlinear solver tolerance and roundoff. We demonstrate exact energy conservation of the method in 3D on the cold two-stream instability.

Yue Zhao, Joshua W. Burby, Andrew Christlieb et al.

We introduce a data-driven approach to learn a generalized kinetic collision operator directly from molecular dynamics. Unlike the conventional (e.g., Landau) models, the present operator takes an anisotropic form that accounts for a second energy transfer arising from the collective interactions between the pair of collision particles and the environment. Numerical results show that preserving the broadly overlooked anisotropic nature of the collision energy transfer is crucial for predicting the plasma kinetics with non-negligible correlations, where the Landau model shows limitations.

Juntao Huang, Yingda Cheng, Andrew J. Christlieb et al.

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

Juntao Huang, Yingda Cheng, Andrew J. Christlieb et al.

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

Juntao Huang, Yingda Cheng, Andrew J. Christlieb et al.

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.

Sami Merhi, Ruochuan Zhang, Mark A. Iwen et al.

In this paper we consider Sparse Fourier Transform (SFT) algorithms for approximately computing the best $s$-term approximation of the Discrete Fourier Transform (DFT) $\mathbf{\hat{f}} \in \mathbb{C}^N$ of any given input vector $\mathbf{f} \in \mathbb{C}^N$ in just $\left( s \log N\right)^{\mathcal{O}(1)}$-time using only a similarly small number of entries of $\mathbf{f}$. In particular, we present a deterministic SFT algorithm which is guaranteed to always recover a near best $s$-term approximation of the DFT of any given input vector $\mathbf{f} \in \mathbb{C}^N$ in $\mathcal{O} \left( s^2 \log ^{\frac{11}{2}} (N) \right)$-time. Unlike previous deterministic results of this kind, our deterministic result holds for both arbitrary vectors $\mathbf{f} \in \mathbb{C}^N$ and vector lengths $N$. In addition to these deterministic SFT results, we also develop several new publicly available randomized SFT implementations for approximately computing $\mathbf{\hat{f}}$ from $\mathbf{f}$ using the same general techniques. The best of these new implementations is shown to outperform existing discrete sparse Fourier transform methods with respect to both runtime and noise robustness for large vector lengths $N$.