Daniel Appelö

Spencer Lee, Daniel Appelö

Fast and accurate classical simulation of quantum systems is a central challenge in the design and control of quantum computers, but the highly oscillatory dynamics of these systems severely limit the efficiency of standard numerical methods. To address this, we adapt Filon quadrature for oscillatory integrals into two numerical methods, called Filon and Controlled Filon, for solving linear systems of ODEs with highly oscillatory solutions. We tailor both methods for efficient implementation in controlled quantum systems, and the Controlled Filon method additionally accounts for the oscillatory structure of the control pulses. We show by numerical experiments that these methods significantly reduce the computational cost of accurately simulating systems of superconducting transmon qubits by decreasing the number of timesteps needed to reach a given level of precision, with only a modest increase in the cost per timestep. For a realistic simulation of the dynamics of a CNOT gate, the Controlled Filon method is the most efficient method tested at every target accuracy, outperforming the best Hermite method by up to 6x and the Hermite method of the same order by up to 500x.

Daniel Appelo, Thomas Hagstrom, Arturo Vargas

Arbitrary order dissipative and conservative Hermite methods for the scalar wave equation achieving $\mathcal{O}(2m)$ orders of accuracy using $\mathcal{O}(m^d)$ degrees of freedom per node in $d$ dimensions are presented. Stability and error analyses as well as implementation strategies for accelerators are also given.

Lu Zhang, Thomas Hagstrom, Daniel Appelo

An energy-based discontinuous Galerkin method for the advective wave equation is proposed and analyzed. Energy-conserving or energy-dissipating methods follow from simple, mesh-independent choices of the inter-element fluxes, and both subsonic and supersonic advection is allowed. Error estimates in the energy norm are established, and numerical experiments on structured grids display optimal convergence in the $L^2$ norm for upwind fluxes. The method generalizes earlier work on energy-based discontinuous Galerkin methods for second order wave equations which was restricted to energy forms written as a simple sum of kinetic and potential energy.

Spencer Lee, Daniel Appelo

This work introduces the High-Order Hermite Optimization (HOHO) method, an open-loop discrete adjoint method for quantum optimal control. Our method is the first of its kind to efficiently compute exact (discrete) gradients when using continuous, parameterized control pulses while solving the forward equations (e.g. Schrodinger's equation or the Linblad master equation) with an arbitrarily high-order Hermite Runge-Kutta method. The HOHO method is implemented in QuantumGateDesign$.$jl (https://github.com/leespen1/QuantumGateDesign.jl), an open-source software package for the Julia programming language, which we use to perform numerical experiments comparing the method to Juqbox$.$jl (https://github.com/LLNL/Juqbox.jl). For realistic model problems we observe speedups up to 775x.

Daniel Appelö, Siyang Wang

We consider wave propagation in a coupled fluid-solid region, separated by a static but possibly curved interface. The wave propagation is modeled by the acoustic wave equation in terms of a velocity potential in the fluid, and the elastic wave equation for the displacement in the solid. At the fluid solid interface, we impose suitable interface conditions to couple the two equations. We use a recently developed, energy based discontinuous Galerkin method to discretize the governing equations in space. Both energy conserving and upwind numerical fluxes are derived to impose the interface conditions. The highlights of the developed scheme include provable energy stability and high order accuracy. We present numerical experiments to illustrate the accuracy property and robustness of the developed scheme.

Tracy Babb, Per-Gunnar Martinsson, Daniel Appelo

A high-order convergent numerical method for solving linear and non-linear parabolic PDEs is presented. The time-stepping is done via an explicit, singly diagonally implicit Runge-Kutta (ESDIRK) method of order 4 or 5, and for the implicit solve, we use the recently developed "Hierarchial Poincare-Steklov (HPS)" method. The HPS method combines a multidomain spectral collocation discretization technique (a "patching method") with a nested-dissection type direct solver. In the context under consideration, the elliptic solve required in each time-step involves the same coefficient matrix, which makes the use of a direct solver particularly effective. The manuscript describes the methodology and presents numerical experiments.

Adeline Kornelus, Daniel Appelö

A new class of Hermite methods for solving nonlinear conservation laws is presented. While preserving the high order spatial accuracy for smooth solutions in the existing Hermite methods, the new methods come with better stability properties. Artificial viscosity in the form of the entropy viscosity method is added to capture shocks.

Mohammad Motamed, Daniel Appelo

We present a Bayesian framework based on a new exponential likelihood function driven by the quadratic Wasserstien metric. Compared to conventional Bayesian models based on Gaussian likelihood functions driven by the least-squares norm ($L_2$ norm), the new framework features several advantages. First, the new framework does not rely on the likelihood of the measurement noise and hence can treat complicated noise structures such as combined additive and multiplicative noise. Secondly, unlike the normal likelihood function, the Wasserstein-based exponential likelihood function does not usually generate multiple local extrema. As a result, the new framework features better convergence to correct posteriors when a Markov Chain Monte Carlo sampling algorithm is employed. Thirdly, in the particular case of signal processing problems, while a normal likelihood function measures only the amplitude differences between the observed and simulated signals, the new likelihood function can capture both the amplitude and the phase differences. We apply the new framework to a class of signal processing problems, that is, the inverse uncertainty quantification of waveforms, and demonstrate its advantages compared to Bayesian models with normal likelihood functions.

Peter DelMastro, Daniel Appelö, Yingda Cheng

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.

Adeline Kornelus, Daniel Appelö

In this work, we outline the entropy viscosity method and discuss how the choice of scaling influences the size of viscosity for a simple shock problem. We present examples to illustrate the performance of the entropy viscosity method under two distinct scalings.

Amit Rotem, Olof Runborg, Daniel Appelö

We consider the numerical solution of the wave equation in materials with rapidly varying coefficients, and time harmonic sources. For these problems, direct discretization is prohibitively costly, and instead multiscale methods are used. There are several multiscale methods that directly discretize in the frequency domain. In this work we instead start in the time-domain and combine a finite difference Heterogeneous Multiscale Method (HMM) for the wave equation with the WaveHoltz method. Each WaveHoltz iteration marches the wave equation towards the time-periodic Helmholtz solution. The advantages of the WaveHoltz method relative to traditional Helmholtz solvers carry over directly to the multiscale problems considered here. Since, in addition, the time-domain solver does not artificially impose boundary conditions on the micro-scale problems, no boundary errors from the micro-scale problems are present in the homogenized frequency domain solution.

Måns I. Andersson, Daniel Appelö

In this paper we explore a fast Poisson solver for problems with a solution that is known to be low-rank. We use an adaptive and warm started cross approximation called Cross-DEIM that iterates between index selection and and cross approximation to generate a low-rank solution. This paper focuses on leveraging a modern machine learning framework, PyTorch, as a general purpose array language to implement low-rank solvers based on Cross-DEIM. PyTorch enables native access to GPUs and accelerators but with a user-friendly high-level interface. We investigate statistical leverage scores for the index selection for the cross approximation due to the cost associated with the pivoted algorithms used with the discrete empirical interpolation methods (DEIM and QDEIM) which are historically preferred. The cross approximation is naturally paired with a Discrete Sine Transform (DST) Poisson solver. This allows the Fast Fourier Transform (FFT) to be evaluated in batches along dimensions independently without any global transpose even in higher dimensions. We present performance results running on a A100 GPU and AMD EPYC CPU demonstrating the usefulness of the approach that enables problems sizes that previously were not feasible.

Daniel Appelo, William D. Henshaw, Zhichao Peng

We present a numerical study of eigenvector deflation as a means of accelerating the WaveHoltz method for solving the Helmholtz equation. For energy-conserving (Dirichlet or Neumann) boundary conditions the WaveHoltz fixed-point iteration converges slowly at high frequency, requiring approximately $\mathcal{O}(ω^{2d})$ iterations in $d$ dimensions. We show that deflating the eigenvectors whose eigenvalues lie nearest the driving frequency substantially reduces iteration counts, and we examine two ways of incorporating the eigenvectors: direct eigenvector deflation (DEVD), in which the forcing and iterate are projected against the deflation set, and augmented-Krylov eigenvector deflation (AUKED) using deflated conjugate gradient (DCG), augmented GMRES (AGMRES), and augmented (recycled) BICGSTAB (ABICGSTAB). The required eigenpairs can be computed efficiently with the EigenWave approach, and we demonstrate, in two dimensions, that when the number of deflation vectors grows quadratically with $ω$ the asymptotic convergence rate remains essentially constant. Because the eigenvectors on structured grids are naturally represented as matrices, we further apply SVD-based compression to reduce their storage. Numerical experiments on single curvilinear grids discretized with summation-by-parts operators, and on overset grids illustrate the robustness and efficiency of the approach, with the deflated solver breaking even against the undeflated solver after as few as two right-hand sides, when accounting for the cost of precomputing the eigenvectors.

Jiuhua Hu, Daniel Appelo, Yingda Cheng

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.