Search

Motivated by applications in unmanned aerial based ground penetratingradar for detecting buried landmines, we consider the problem of imagingsmall point like scatterers situated in a lossy medium below a randomrough surface. Both the random rough surface and the absorption in thelossy medium significantly impede the target detection and imagingprocess. Using principal component analysis we effectively remove thereflection from the air-soil interface. We then use a modification ofthe classical synthetic aperture radar imaging functional to image thetargets. This imaging method introduces a user-defined parameter, δ,which scales the resolution by √δ allowing for target localization withsub wavelength accuracy. Numerical results in two dimensionsillustrateWe study imaging methods for identifying point targets in alossy medium below a random rough surface. the robustness of theapproach for imaging multiple targets. However, the depth at whichtargets are detectable is limited due to the absorption in the lossymedium.

We introduce a dispersive point target model based on scattering by aparticle in the far-field. The synthetic aperture imaging problem is thenexpanded to identify these targets and recover their locations and frequencydependent reflectivities. We show that Kirchhoff migration (KM) is able toidentify dispersive point targets in an imaging region. However, KM predictstarget locations that are shifted in range from their true locations. We derivean estimate for this range shift for a single target. We also show that becauseof this range shift we cannot recover the complex-valued frequency dependentreflectivity, but we can recover its absolute value and hence the radarcross-section (RCS) of the target. Simulation results show that we can detect,recover the approximate location, and recover the RCS for dispersive pointtargets thereby opening opportunities to classifying important differencesbetween multiple targets such as their sizes or material compositions.

We present four frequently used finite difference methods and establish the error bounds for the discretization of the Dirac equation in the massless and nonrelativistic regime, involving a small dimensionless parameter 0 < ε ≪ 1 inversely proportional to the speed of light. In the massless and nonrelativistic regime, the solution exhibits rapid motion in space and is highly oscillatory in time. Specifically, the wavelength of the propagating waves in time is at O(ε), while in space, it is at O(1) with the wave speed at O(ε− 1). We adopt one leap-frog, two semi-implicit, and one conservative Crank-Nicolson finite difference methods to numerically discretize the Dirac equation in one dimension and establish rigorously the error estimates which depend explicitly on the time step τ, mesh size h, and the small parameter ε. The error bounds indicate that, to obtain the “correct” numerical solution in the massless and nonrelativistic regime, i.e., 0 < ε ≪ 1, all these finite difference methods share the same ε-scalability as time step τ = O(ε3/2) and mesh size h = O(ε1/2). A large number of numerical results are reported to verify the error estimates.

The goal of the present article is to understand the impact of numerical schemes for the reconstruction of data at cell faces in finite-volume methods, and to assess their interaction with the quadrature rule used to compute the average over the cell volume. Here, third-, fifth- and seventh-order WENO-Z schemes are investigated. On a problem with a smooth solution, the theoretical order of convergence rate for each method is retrieved, and changing the order of the reconstruction at cell faces does not impact the results, whereas for a shock-driven problem all the methods collapse to first-order. Study of the decay of compressible homogeneous isotropic turbulence reveals that using a high-order quadrature rule to compute the average over a finite-volume cell does not improve the spectral accuracy and that all methods present a second-order convergence rate. However the choice of the numerical method to reconstruct data at cell faces is found to be critical to correctly capture turbulent spectra. In the context of simulations with finite-volume methods of practical flows encountered in engineering applications, it becomes apparent that an efficient strategy is to perform the average integration with a low-order quadrature rule on a fine mesh resolution, whereas high-order schemes should be used to reconstruct data at cell faces.

Quantum computing is a promising technology that harnesses the peculiaritiesof quantum mechanics to deliver computational speedups for some problems thatare intractable to solve on a classical computer. Current generation noisyintermediate-scale quantum (NISQ) computers are severely limited in terms ofchip size and error rates. Shallow quantum circuits with uncomplicatedtopologies are essential for successful applications in the NISQ era. Based onmatrix analysis, we derive localized circuit transformations to efficientlycompress quantum circuits for simulation of certain spin Hamiltonians known asfree fermions. The depth of the compressed circuits is independent ofsimulation time and grows linearly with the number of spins. The proposednumerical circuit compression algorithm behaves backward stable and scalescubically in the number of spins enabling circuit synthesis beyond$\mathcal{O}(10^3)$ spins. The resulting quantum circuits have a simplenearest-neighbor topology, which makes them ideally suited for NISQ devices.

Exponential integrators are a well-known class of time integration methodsthat have been the subject of many studies and developments in the past twodecades. Surprisingly, there have been limited efforts to analyze theirstability and efficiency on non-diffusive equations to date. In this paper weapply linear stability analysis to showcase the poor stability properties ofexponential integrators on non-diffusive problems. We then propose a simplerepartitioning approach that stabilizes the integrators and enables theefficient solution of stiff, non-diffusive equations. To validate theeffectiveness of our approach, we perform several numerical experiments thatcompare partitioned exponential integrators to unmodified ones. We also comparerepartitioning to the well-known approach of adding hyperviscosity to theequation right-hand-side. Overall, we find that the repartitioning restoresconvergence at large timesteps and, unlike hyperviscosity, it does not requirethe use of high-order spatial derivatives.

We introduce a set of novel multiscale basis transforms for signals on graphs that utilize their “dual” domains by incorporating the “natural” distances between graph Laplacian eigenvectors, rather than simply using the eigenvalue ordering. These basis dictionaries can be seen as generalizations of the classical Shannon wavelet packet dictionary to arbitrary graphs, and do not rely on the frequency interpretation of Laplacian eigenvalues. We describe the algorithms (involving either vector rotations or orthogonalizations) to construct these basis dictionaries, use them to efficiently approximate graph signals through the best basis search, and demonstrate the strengths of these basis dictionaries for graph signals measured on sunflower graphs and street networks.

We present a computational framework for solving the equations of inviscidgas dynamics using structured grids with embedded geometries. The novelty ofthe proposed approach is the use of high-order discontinuous Galerkin (dG)schemes and a shock-capturing Finite Volume (FV) scheme coupled via an $hp$adaptive mesh refinement ($hp$-AMR) strategy that offers high-order accurateresolution of the embedded geometries. The $hp$-AMR strategy is based on amulti-level block-structured domain partition in which each level isrepresented by block-structured Cartesian grids and the embedded geometry isrepresented implicitly by a level set function. The intersection of theembedded geometry with the grids produces the implicitly-defined mesh thatconsists of a collection of regular rectangular cells plus a relatively smallnumber of irregular curved elements in the vicinity of the embedded boundaries.High-order quadrature rules for implicitly-defined domains enable high-orderaccuracy resolution of the curved elements with a cell-merging strategy toaddress the small-cell problem. The $hp$-AMR algorithm treats the system with asecond-order finite volume scheme at the finest level to dynamically track theevolution of solution discontinuities while using dG schemes at coarser levelsto provide high-order accuracy in smooth regions of the flow. On the dG levels,the methodology supports different orders of basis functions on differentlevels. The space-discretized governing equations are then advanced explicitlyin time using high-order Runge-Kutta algorithms. Numerical tests are presentedfor two-dimensional and three-dimensional problems involving an ideal gas. Theresults are compared with both analytical solutions and experimentalobservations and demonstrate that the framework provides high-order accuracyfor smooth flows and accurately captures solution discontinuities.

We compare the long-time error bounds and spatial resolution of finitedifference methods with different spatial discretizations for the Diracequation with small electromagnetic potentials characterized by $\varepsilon\in (0, 1]$ a dimensionless parameter. We begin with the simple and widely usedfinite difference time domain (FDTD) methods, and establish rigorous errorbounds of them, which are valid up to the time at $O(1/\varepsilon)$. In theerror estimates, we pay particular attention to how the errors dependexplicitly on the mesh size $h$ and time step $\tau$ as well as the smallparameter $\varepsilon$. Based on the results, in order to obtain "correct"numerical solutions up to the time at $O(1/\varepsilon)$, the$\varepsilon$-scalability (or meshing strategy requirement) of the FDTD methodsshould be taken as $h = O(\varepsilon^{1/2})$ and $\tau =O(\varepsilon^{1/2})$. To improve the spatial resolution capacity, we apply theFourier spectral method to discretize the Dirac equation in space. Error boundsof the resulting finite difference Fourier pseudospectral (FDFP) methods showthat they exhibit uniform spatial errors in the long-time regime, which areoptimal in space as suggested by the Shannon's sampling theorem. Extensivenumerical results are reported to confirm the error bounds and demonstrate thatthey are sharp.

The Hartree-Fock-Bogoliubov (HFB) theory is the starting point for treating superconducting systems. However, the computational cost for solving large scale HFB equations can be much larger than that of the Hartree-Fock equations, particularly when the Hamiltonian matrix is sparse, and the number of electrons N is relatively small compared to the matrix size Nb. We first provide a concise and relatively self-contained review of the HFB theory for general finite sized quantum systems, with special focus on the treatment of spin symmetries from a linear algebra perspective. We then demonstrate that the pole expansion and selected inversion (PEXSI) method can be particularly well suited for solving large scale HFB equations. For a Hubbard-Type Hamiltonian, the cost of PEXSI is at most (Nb2) for both gapped and gapless systems, which can be significantly faster than the standard cubic scaling diagonalization methods. We show that PEXSI can solve a two-dimensional Hubbard-Hofstadter model with Nb up to 2.88 × 106, and the wall clock time is less than 100 s using 17 280 CPU cores. This enables the simulation of physical systems under experimentally realizable magnetic fields, which cannot be otherwise simulated with smaller systems.

In this paper, we present an approach to deal with the dynamics of the Diracequation with time-dependent electromagnetic potentials using the fourth-ordercompact time-splitting method ($S_\text{4c}$). To this purpose, thetime-ordering technique for time-dependent Hamiltonians is introduced, so thatthe influence of the time-dependence could be limited to certain steps whichare easy to treat. Actually, in the case of the Dirac equation, it turns outthat only those steps involving potentials need to be amended, and the schemeremains efficient, accurate, as well as easy to implement. Numerical examplesin 1D and 2D are given to validate the scheme.

We present a numerical method specifically designed for simulatingthree-dimensional fluid--structure interaction (FSI) problems based on thereference map technique (RMT). The RMT is a fully Eulerian FSI numerical methodthat allows fluids and large-deformation elastic solids to be represented on asingle fixed computational grid. This eliminates the need for meshing complexgeometries typical in other FSI approaches, and greatly simplifies the couplingbetween fluid and solids. We develop the first three-dimensional implementationof the RMT, parallelized using the distributed memory paradigm, to simulateincompressible FSI with neo-Hookean solids. As part of our new method, wedevelop a new field extrapolation scheme that works efficiently in parallel.Through representative examples, we demonstrate the method's accuracy andconvergence, as well as its suitability in investigating many-body and activesystems. The examples include settling of a mixture of heavy and buoyant softellipsoids, lid-driven cavity flow containing a soft sphere, and swimmersactuated via active stress.

The density-equalizing map, a technique developed for cartogram creation, has been widely applied to data visualization but only for 2D applications. In this work, we propose a novel method called the volumetric density-equalizing reference map for computing density-equalizing map for volumetric domains. Given a prescribed density distribution in a volumetric domain in R3, the proposed method continuously deforms the domain, with different volume elements enlarged or shrunk according to the density distribution. With the aid of the proposed method, medical and sociological data can be visualized via deformations of 3D objects. The method can also be applied to adaptive remeshing and shape modeling. Furthermore, by exploiting the time-dependent nature of the proposed method, applications to shape morphing can be easily achieved. Experimental results are presented to demonstrate the effectiveness of the proposed method.

An extension of the two-step staggered time discretization of linear elastodynamics in stress-velocity form to systems involving internal variables subjected to a possibly non-linear dissipative evolution is proposed. The original scheme is thus enhanced by another step for the internal variables which, in general, is implicit, although even this step might be explicit if no spatial gradients of the internal variables are involved. Using an abstract Banach-space formulation, a priori estimates and convergence are proved under a CFL condition. The developed three-step scheme finds applications in various problems of continuum mechanics at small strain. Here, we consider in particular plasticity, viscoelasticity (creep), diffusion in poroelastic media, and damage.

The Fragile Points method (FPM) is an elementarily simple Galerkin meshless method, employing Point-based discontinuous trial and test functions only, without using element-based trial and test functions. In this study, the algorithmic formulations of FPM for linear elasticity are given in detail, by exploring the concepts of point stiffness matrices and numerical flux corrections. Advantages of FPM for simulating the deformations of complex structures, and for simulating complex crack propagations and rupture developments, are also thoroughly discussed. Numerical examples of deformation and stress analyses of benchmark problems, as well as of realistic structures with complex geometries, demonstrate the accuracy, efficiency, and robustness of the proposed FPM. Simulations of crack initiation and propagations are also given in this study, demonstrating the advantages of the present FPM in modeling complex rupture and fracture phenomena. The crack and rupture propagation modeling in FPM is achieved without remeshing or augmenting the trial functions as in standard, extended, or generalized finite element method. The simulation of impact, penetration, and other extreme problems by FPM will be discussed in our future articles.

We examine a method for solving an infinite-dimensional tensor eigenvalueproblem $H x = \lambda x$, where the infinite-dimensional symmetric matrix $H$exhibits a translational invariant structure. We provide a formulation of thistype of problem from a numerical linear algebra point of view and describe howa power method applied to $e^{-Ht}$ is used to obtain an approximation to thedesired eigenvector. This infinite-dimensional eigenvector is represented in acompact way by a translational invariant infinite Tensor Ring (iTR). Low rankapproximation is used to keep the cost of subsequent power iterations boundedwhile preserving the iTR structure of the approximate eigenvector. We show howthe averaged Rayleigh quotient of an iTR eigenvector approximation can beefficiently computed and introduce a projected residual to monitor itsconvergence. In the numerical examples, we illustrate that the norm of thisprojected iTR residual can also be used to automatically modify the time step$t$ to ensure accurate and rapid convergence of the power method.

We propose a flexible power method for computing the leftmost, i.e.,algebraically smallest, eigenvalue of an infinite dimensional tensor eigenvalueproblem, $H x = \lambda x$, where the infinite dimensional symmetric matrix $H$exhibits a translational invariant structure. We assume the smallest eigenvalueof $H$ is simple and apply a power iteration of $e^{-H}$ with the eigenvectorrepresented in a compact way as a translational invariant infinite Tensor Ring(iTR). Hence, the infinite dimensional eigenvector can be represented by afinite number of iTR cores of finite rank. In order to implement this poweriteration, we use a small parameter $t$ so that the infinite matrix-vectoroperation $e^{-Ht}x$ can efficiently be approximated by the Lie productformula, also known as Suzuki--Trotter splitting, and we employ a low rankapproximation through a truncated singular value decomposition on the iTR coresin order to keep the cost of subsequent power iterations bounded. We also usean efficient way for computing the iTR Rayleigh quotient and introduce a finitesize iTR residual which is used to monitor the convergence of the Rayleighquotient and to modify the timestep $t$. In this paper, we discuss 2 differentimplementations of the flexible power algorithm and illustrate the automatictimestep adaption approach for several numerical examples.

This paper presents an adaptive randomized algorithm for computing the butterfly factorization of an m × n matrix with m ≈ n provided that both the matrix and its transpose can be rapidly applied to arbitrary vectors. The resulting factorization is composed of O(log n) sparse factors, each containing O(n) nonzero entries. The factorization can be attained using O(n3/2 log n) computation and O(n log n) memory resources. The proposed algorithm can be implemented in parallel and can apply to matrices with strong or weak admissibility conditions arising from surface integral equation solvers as well as multi-frontal-based finite-difference, finite-element, or finite-volume solvers. A distributed-memory parallel implementation of the algorithm demonstrates excellent scaling behavior.

The fast Fourier transform (FFT) is one of the most successful numerical algorithms of the 20th century and has found numerous applications in many branches of computational science and engineering. The FFT algorithm can be derived from a particular matrix decomposition of the discrete Fourier transform (DFT) matrix. In this paper, we show that the quantum Fourier transform (QFT) can be derived by further decomposing the diagonal factors of the FFT matrix decomposition into products of matrices with Kronecker product structure. We analyze the implication of this Kronecker product structure on the discrete Fourier transform of rank-1 tensors on a classical computer. We also explain why such a structure can take advantage of an important quantum computer feature that enables the QFT algorithm to attain an exponential speedup on a quantum computer over the FFT algorithm on a classical computer. Further, the connection between the matrix decomposition of the DFT matrix and a quantum circuit is made. We also discuss a natural extension of a radix-2 QFT decomposition to a radix-d QFT decomposition. No prior knowledge of quantum computing is required to understand what is presented in this paper. Yet, we believe this paper may help readers to gain some rudimentary understanding of the nature of quantum computing from a matrix computation point of view.

We present four frequently used finite difference methods and establish the error bounds for the discretization of the Dirac equation in the massless and nonrelativistic regime, involving a small dimensionless parameter 0 < ε ≪ 1 inversely proportional to the speed of light. In the massless and nonrelativistic regime, the solution exhibits rapid motion in space and is highly oscillatory in time. Specifically, the wavelength of the propagating waves in time is at O(ε), while in space, it is at O(1) with the wave speed at O(ε− 1). We adopt one leap-frog, two semi-implicit, and one conservative Crank-Nicolson finite difference methods to numerically discretize the Dirac equation in one dimension and establish rigorously the error estimates which depend explicitly on the time step τ, mesh size h, and the small parameter ε. The error bounds indicate that, to obtain the “correct” numerical solution in the massless and nonrelativistic regime, i.e., 0 < ε ≪ 1, all these finite difference methods share the same ε-scalability as time step τ = O(ε3/2) and mesh size h = O(ε1/2). A large number of numerical results are reported to verify the error estimates.

The central importance of large-scale eigenvalue problems in scientific computation necessitates the development of massively parallel algorithms for their solution. Recent advances in dense numerical linear algebra have enabled the routine treatment of eigenvalue problems with dimensions on the order of hundreds of thousands on the world's largest supercomputers. In cases where dense treatments are not feasible, Krylov subspace methods offer an attractive alternative due to the fact that they do not require storage of the problem matrices. However, demonstration of scalability of either of these classes of eigenvalue algorithms on computing architectures capable of expressing massive parallelism is non-trivial due to communication requirements and serial bottlenecks, respectively. In this work, we introduce the SISLICE method: a parallel shift-invert algorithm for the solution of the symmetric self-consistent field (SCF) eigenvalue problem. The SISLICE method drastically reduces the communication requirement of current parallel shift-invert eigenvalue algorithms through various shift selection and migration techniques based on density of states estimation and k-means clustering, respectively. This work demonstrates the robustness and parallel performance of the SISLICE method on a representative set of SCF eigenvalue problems and outlines research directions that will be explored in future work.

Parareal is a widely studied parallel-in-time method that can achievemeaningful speedup on certain problems. However, it is well known that themethod typically performs poorly on non-diffusive equations. This paperanalyzes linear stability and convergence for IMEX Runge-Kutta Parareal methodson non-diffusive equations. By combining standard linear stability analysiswith a simple convergence analysis, we find that certain Pararealconfigurations can achieve parallel speedup on non-diffusive equations. Thesestable configurations all posses low iteration counts, large block sizes, and alarge number of processors. Numerical examples using the nonlinear Schrodingerequation demonstrate the analytical conclusions.

We present a general simulation approach for fluid-solid interactions based on the fully Eulerian reference map technique. The approach permits the modelling of one or more finitely deformable continuum solid bodies interacting with a fluid and with each other. A key advantage of this approach is its ease of use, as the solid and fluid are discretized on the same fixed grid, which greatly simplifies the coupling between the phases. We use the method to study a number of illustrative examples involving an incompressible Navier-Stokes fluid interacting with multiple neo-Hookean solids. Our method has several useful features including the ability to model solids with sharp corners and the ability to model actuated solids. The latter permits the simulation of active media such as swimmers, which we demonstrate. The method is validated favourably in the flag-flapping geometry, for which a number of experimental, numerical and analytical studies have been performed. We extend the flapping analysis beyond the thin-flag limit, revealing an additional destabilization mechanism to induce flapping.

This article presents a low-rank decomposition algorithm based on subsampling of matrix entries. The proposed algorithm first computes rank-revealing decompositions of submatrices with a blocked adaptive cross approximation (BACA) algorithm, and then applies a hierarchical merge operation via truncated singular value decompositions (H-BACA). The proposed algorithm significantly improves the convergence of the baseline ACA algorithm and achieves reduced computational complexity compared to the traditional decompositions such as rank-revealing QR. Numerical results demonstrate the efficiency, accuracy, and parallel scalability of the proposed algorithm.

We utilize generalized moving least squares (GMLS) to develop meshfreetechniques for discretizing hydrodynamic flow problems on manifolds. We useexterior calculus to formulate incompressible hydrodynamic equations in theStokesian regime and handle the divergence-free constraints via a generalizedvector potential. This provides less coordinate-centric descriptions andenables the development of efficient numerical methods and splitting schemesfor the fourth-order governing equations in terms of a system of second-orderelliptic operators. Using a Hodge decomposition, we develop methods formanifolds having spherical topology. We show the methods exhibit high-orderconvergence rates for solving hydrodynamic flows on curved surfaces. Themethods also provide general high-order approximations for the metric,curvature, and other geometric quantities of the manifold and associatedexterior calculus operators. The approaches also can be utilized to develophigh-order solvers for other scalar-valued and vector-valued problems onmanifolds.

Hyperbolic systems under nonconservative form arise in numerous applicationsmodeling physical processes, for example from the relaxation of more generalequations (e.g. with dissipative terms). This paper reviews an existing classof augmented Roe schemes and discusses their application to linearnonconservative hyperbolic systems with source terms. We extend existingaugmented methods by redefining them within a common framework which uses ageometric reinterpretation of source terms. This results in intrinsicallywell-balanced numerical discretizations. We discuss two equivalentformulations: (1) a nonconservative approach and (2) a conservativereformulation of the problem. The equilibrium properties of the schemes areexamined and the conditions for the preservation of the well-balanced propertyare provided. Transient and steady state test cases for linear acoustics andhyperbolic heat equations are presented. A complete set of benchmark problemswith analytical solution, including transient and steady situations withdiscontinuities in the medium properties, are presented and used to assess theequilibrium properties of the schemes. It is shown that the proposed schemessatisfy the expected equilibrium and convergence properties.

Mean field games (MFG) and mean field control (MFC) are critical classes of multiagent models for the efficient analysis of massive populations of interacting agents. Their areas of application span topics in economics, finance, game theory, industrial engineering, crowd motion, and more. In this paper, we provide a flexible machine learning framework for the numerical solution of potential MFG and MFC models. State-of-the-art numerical methods for solving such problems utilize spatial discretization that leads to a curse of dimensionality. We approximately solve high-dimensional problems by combining Lagrangian and Eulerian viewpoints and leveraging recent advances from machine learning. More precisely, we work with a Lagrangian formulation of the problem and enforce the underlying Hamilton-Jacobi-Bellman (HJB) equation that is derived from the Eulerian formulation. Finally, a tailored neural network parameterization of the MFG/MFC solution helps us avoid any spatial discretization. Our numerical results include the approximate solution of 100-dimensional instances of optimal transport and crowd motion problems on a standard work station and a validation using a Eulerian solver in two dimensions. These results open the door to much-anticipated applications of MFG and MFC models that are beyond reach with existing numerical methods.

The modeling of atmospheric processes in the context of weather and climate simulations is an important and computationally expensive challenge. The temporal integration of the underlying PDEs requires a very large number of time steps, even when the terms accounting for the propagation of fast atmospheric waves are treated implicitly. Therefore, the use of parallel-in-time integration schemes to reduce the time-to-solution is of increasing interest, particularly in the numerical weather forecasting field. We present a multi-level parallel-in-time integration method combining the Parallel Full Approximation Scheme in Space and Time (PFASST) with a spatial discretization based on Spherical Harmonics (SH). The iterative algorithm computes multiple time steps concurrently by interweaving parallel high-order fine corrections and serial corrections performed on a coarsened problem. To do that, we design a methodology relying on the spectral basis of the SH to coarsen and interpolate the problem in space. The methods are evaluated on the shallow-water equations on the sphere using a set of tests commonly used in the atmospheric flow community. We assess the convergence of PFASST-SH upon refinement in time. We also investigate the impact of the coarsening strategy on the accuracy of the scheme, and specifically on its ability to capture the high-frequency modes accumulating in the solution. Finally, we study the computational cost of PFASST-SH to demonstrate that our scheme resolves the main features of the solution multiple times faster than the serial schemes.

The modeling of atmospheric processes in the context of weather and climate simulations is an important and computationally expensive challenge. The temporal integration of the underlying PDEs requires a very large number of time steps, even when the terms accounting for the propagation of fast atmospheric waves are treated implicitly. Therefore, the use of parallel-in-time integration schemes to reduce the time-to-solution is of increasing interest, particularly in the numerical weather forecasting field. We present a multi-level parallel-in-time integration method combining the Parallel Full Approximation Scheme in Space and Time (PFASST) with a spatial discretization based on Spherical Harmonics (SH). The iterative algorithm computes multiple time steps concurrently by interweaving parallel high-order fine corrections and serial corrections performed on a coarsened problem. To do that, we design a methodology relying on the spectral basis of the SH to coarsen and interpolate the problem in space. The methods are evaluated on the shallow-water equations on the sphere using a set of tests commonly used in the atmospheric flow community. We assess the convergence of PFASST-SH upon refinement in time. We also investigate the impact of the coarsening strategy on the accuracy of the scheme, and specifically on its ability to capture the high-frequency modes accumulating in the solution. Finally, we study the computational cost of PFASST-SH to demonstrate that our scheme resolves the main features of the solution multiple times faster than the serial schemes.

Starting from measured data, we develop a method to compute the finestructure of the spectrum of the Koopman operator with rigorous convergenceguarantees. The method is based on the observation that, in themeasure-preserving ergodic setting, the moments of the spectral measureassociated to a given observable are computable from a single trajectory ofthis observable. Having finitely many moments available, we use the classicalChristoffel-Darboux kernel to separate the atomic and absolutely continuousparts of the spectrum, supported by convergence guarantees as the number ofmoments tends to infinity. In addition, we propose a technique to detect thesingular continuous part of the spectrum as well as two methods to approximatethe spectral measure with guaranteed convergence in the weak topology,irrespective of whether the singular continuous part is present or not. Theproposed method is simple to implement and readily applicable to large-scalesystems since the computational complexity is dominated by inverting an$N\times N$ Hermitian positive-definite Toeplitz matrix, where $N$ is thenumber of moments, for which efficient and numerically stable algorithms exist;in particular, the complexity of the approach is independent of the dimensionof the underlying state-space. We also show how to compute, from measured data,the spectral projection on a given segment of the unit circle, allowing us toobtain a finite-dimensional approximation of the operator that explicitly takesinto account the point and continuous parts of the spectrum. Finally, wedescribe a relationship between the proposed method and the so-called HankelDynamic Mode Decomposition, providing new insights into the behavior of theeigenvalues of the Hankel DMD operator. A number of numerical examplesillustrate the approach, including a study of the spectrum of the lid-driventwo-dimensional cavity flow.

We consider the problem of imaging sparse scenes from a few noisy data usingan $l_1$-minimization approach. This problem can be cast as a linear system ofthe form $A \, \rho =b$, where $A$ is an $N\times K$ measurement matrix. Weassume that the dimension of the unknown sparse vector $\rho \in{\mathbb{C}}^K$ is much larger than the dimension of the data vector $b \in{\mathbb{C}}^N$, i.e, $K \gg N$. We provide a theoretical framework that allowsus to examine under what conditions the $\ell_1$-minimization problem admits asolution that is close to the exact one in the presence of noise. Our analysisshows that $l_1$-minimization is not robust for imaging with noisy data whenhigh resolution is required. To improve the performance of $l_1$-minimizationwe propose to solve instead the augmented linear system $ [A \, | \, C] \rho=b$, where the $N \times \Sigma$ matrix $C$ is a noise collector. It isconstructed so as its column vectors provide a frame on which the noise of thedata, a vector of dimension $N$, can be well approximated. Theoretically, thedimension $\Sigma$ of the noise collector should be $e^N$ which would make itsuse not practical. However, our numerical results illustrate that robustresults in the presence of noise can be obtained with a large enough number ofcolumns $\Sigma \approx 10 K$.

In this paper we analyze a finite element method applied to a continuous downscaling data assimilation algorithm for the numerical approximation of the two- and three-dimensional Navier-Stokes equations corresponding to given measurements on a coarse spatial scale. For representing the coarse mesh measurements we consider different types of interpolation operators including a Lagrange interpolant. We obtain uniform-in-time estimates for the error between a finite element approximation and the reference solution corresponding to the coarse mesh measurements. We consider both the case of a plain Galerkin method and a Galerkin method with grad-div stabilization. For the stabilized method we prove error bounds in which the constants do not depend on inverse powers of the viscosity. Some numerical experiments illustrate the theoretical results.

We consider imaging the reflectivity of scatterers from intensity-only datarecorded by a single moving transducer that both emits and receives signals,forming a synthetic aperture. By exploiting frequency illumination diversity,we obtain multiple intensity measurements at each location, from which wedetermine field cross-correlations using an appropriate phase controlledillumination strategy and the inner product polarization identity. The fieldcross-correlations obtained this way do not, however, provide all the missingphase information because they are determined up to a phase that depends on thereceiver's location. The main result of this paper is an algorithm with whichwe recover the field cross-correlations up to a single phase that is common toall the data measured over the synthetic aperture, so all the data aresynchronized. Thus, we can image coherently with data over all frequencies andmeasurement locations as if full phase information was recorded.

Pole-swapping algorithms, which are generalizations of the QZ algorithm for the generalized eigenvalue problem, are studied. A new modular (and therefore more flexible) convergence theory that applies to all pole-swapping algorithms is developed. A key component of all such algorithms is a procedure that swaps two adjacent eigenvalues in a triangular pencil. An improved swapping routine is developed, and its superiority over existing methods is demonstrated by a backward error analysis and numerical tests. The modularity of the new convergence theory and the generality of the pole-swapping approach shed new light on bi-directional chasing algorithms, optimally packed shifts, and bulge pencils, and allow the design of novel algorithms.

We present a general simulation approach for fluid-solid interactions based on the fully Eulerian reference map technique. The approach permits the modelling of one or more finitely deformable continuum solid bodies interacting with a fluid and with each other. A key advantage of this approach is its ease of use, as the solid and fluid are discretized on the same fixed grid, which greatly simplifies the coupling between the phases. We use the method to study a number of illustrative examples involving an incompressible Navier-Stokes fluid interacting with multiple neo-Hookean solids. Our method has several useful features including the ability to model solids with sharp corners and the ability to model actuated solids. The latter permits the simulation of active media such as swimmers, which we demonstrate. The method is validated favourably in the flag-flapping geometry, for which a number of experimental, numerical and analytical studies have been performed. We extend the flapping analysis beyond the thin-flag limit, revealing an additional destabilization mechanism to induce flapping.

Pole-swapping algorithms, which are generalizations of the QZ algorithm for the generalized eigenvalue problem, are studied. A new modular (and therefore more flexible) convergence theory that applies to all pole-swapping algorithms is developed. A key component of all such algorithms is a procedure that swaps two adjacent eigenvalues in a triangular pencil. An improved swapping routine is developed, and its superiority over existing methods is demonstrated by a backward error analysis and numerical tests. The modularity of the new convergence theory and the generality of the pole-swapping approach shed new light on bi-directional chasing algorithms, optimally packed shifts, and bulge pencils, and allow the design of novel algorithms.

A unified construction of the Hm-nonconforming virtual elements of any order k is developed on any shape of polytope in Rn with constraints m ≤ n and k ≥ m. As a vital tool in the construction, a generalized Green's identity for Hm inner product is derived. The Hm-nonconforming virtual element methods are then used to approximate solutions of the m-harmonic equation. After establishing a bound on the jump related to the weak continuity, the optimal error estimate of the canonical interpolation, and the norm equivalence of the stabilization term, the optimal error estimates are derived for the Hm-nonconforming virtual element methods.

We introduce a general framework for approximating parabolic StochasticPartial Differential Equations (SPDEs) based on fluctuation-dissipationbalance. Using this approach we formulate Stochastic Discontinuous GalerkinMethods (SDGM). We show how methods with linear-time computational complexitycan be developed for handling domains with general geometry and generatingstochastic terms handling both Dirichlet and Neumann boundary conditions. Wedemonstrate our approach on example systems and contrast with alternativeapproaches using direct stochastic discretizations based on random fluxes. Weshow how our Fluctuation-Dissipation Discretizations (FDD) framework allows forcompensating for differences in dissipative properties of discrete numericaloperators relative to their continuum counter-parts. This allows us to handlegeneral heterogeneous discretizations capturing accurately statisticalrelations. Our FDD framework provides a general approach for formulating SDGMdiscretizations and other numerical methods for robust approximation ofstochastic differential equations.

In recent literature, a general two step procedure has been formulated forsolving the problem of phase retrieval. First, a spectral technique is used toobtain a constant-error initial estimate, following which, the estimate isrefined to arbitrary precision by first-order optimization of a non-convex lossfunction. Numerical experiments, however, seem to suggest that simply runningthe iterative schemes from a random initialization may also lead toconvergence, albeit at the cost of slightly higher sample complexity. In thispaper, we prove that, in fact, constant step size online stochastic gradientdescent (SGD) converges from arbitrary initializations for the non-smooth,non-convex amplitude squared loss objective. In this setting, online SGD isalso equivalent to the randomized Kaczmarz algorithm from numerical analysis.Our analysis can easily be generalized to other single index models. It alsomakes use of new ideas from stochastic process theory, including the notion ofa summary state space, which we believe will be of use for the broader field ofnon-convex optimization.

This paper presents a general theory and isogeometric finite element implementation for studying mass conserving phase transitions on deforming surfaces. The mathematical problem is governed by two coupled fourth-order nonlinear partial differential equations (PDEs) that live on an evolving two-dimensional manifold. For the phase transitions, the PDE is the Cahn–Hilliard equation for curved surfaces, which can be derived from surface mass balance in the framework of irreversible thermodynamics. For the surface deformation, the PDE is the (vector-valued) Kirchhoff–Love thin shell equation. Both PDEs can be efficiently discretized using C 1 -continuous interpolations without derivative degrees-of-freedom (dofs). Structured NURBS and unstructured spline spaces with pointwise C 1 -continuity are utilized for these interpolations. The resulting finite element formulation is discretized in time by the generalized-α scheme with adaptive time-stepping, and it is fully linearized within a monolithic Newton–Raphson approach. A curvilinear surface parameterization is used throughout the formulation to admit general surface shapes and deformations. The behavior of the coupled system is illustrated by several numerical examples exhibiting phase transitions on deforming spheres, tori and double-tori.

First-principles Markov Chain Monte Carlo sampling is used to investigate uncertainty quantification and uncertainty propagation in parameters describing hydrogen kinetics. Specifically, we sample the posterior distribution for thirty-one parameters focusing on the H2O2 and HO2 reactions resulting from conditioning on ninety-one experiments. Established literature values are used for the remaining parameters in the mechanism as well as other thermodynamic and transport data needed to specify fluid properties. The samples are computed using an affine invariant sampler starting with broad, noninformative priors. Autocorrelation analysis shows that O(1M) samples are sufficient to obtain a reasonable sampling of the posterior. The resulting distribution identifies strong positive and negative correlations and several non-Gaussian characteristics. Using samples drawn from the posterior, we investigate the impact of parameter uncertainty on the prediction of two more complex flames: a 2D premixed flame kernel and the ignition of a hydrogen jet issuing into a heated chamber. The former represents a combustion regime similar to the target experiments used to calibrate the mechanism and the latter represents a different combustion regime. For the premixed flame, the net amount of product after a given time interval has a standard deviation of less than 2% whereas the standard deviation of the ignition time for the jet is more than 10%. The samples used for these studies are posted online. These results indicate the degree to which parameters consistent with the target experiments constrain predicted behavior in different combustion regimes. This process provides a framework for both identifying reactions for further study from candidate mechanisms as well as combining uncertainty quantification and propagation to, ultimately, tie uncertainty in laboratory flame experiments to uncertainty in end-use numerical predictions of more complicated scenarios.

This paper presents a general theory and isogeometric finite element implementation for studying mass conserving phase transitions on deforming surfaces. The mathematical problem is governed by two coupled fourth-order nonlinear partial differential equations (PDEs) that live on an evolving two-dimensional manifold. For the phase transitions, the PDE is the Cahn–Hilliard equation for curved surfaces, which can be derived from surface mass balance in the framework of irreversible thermodynamics. For the surface deformation, the PDE is the (vector-valued) Kirchhoff–Love thin shell equation. Both PDEs can be efficiently discretized using C 1 -continuous interpolations without derivative degrees-of-freedom (dofs). Structured NURBS and unstructured spline spaces with pointwise C 1 -continuity are utilized for these interpolations. The resulting finite element formulation is discretized in time by the generalized-α scheme with adaptive time-stepping, and it is fully linearized within a monolithic Newton–Raphson approach. A curvilinear surface parameterization is used throughout the formulation to admit general surface shapes and deformations. The behavior of the coupled system is illustrated by several numerical examples exhibiting phase transitions on deforming spheres, tori and double-tori.

First-principles Markov Chain Monte Carlo sampling is used to investigate uncertainty quantification and uncertainty propagation in parameters describing hydrogen kinetics. Specifically, we sample the posterior distribution for thirty-one parameters focusing on the H2O2 and HO2 reactions resulting from conditioning on ninety-one experiments. Established literature values are used for the remaining parameters in the mechanism as well as other thermodynamic and transport data needed to specify fluid properties. The samples are computed using an affine invariant sampler starting with broad, noninformative priors. Autocorrelation analysis shows that O(1M) samples are sufficient to obtain a reasonable sampling of the posterior. The resulting distribution identifies strong positive and negative correlations and several non-Gaussian characteristics. Using samples drawn from the posterior, we investigate the impact of parameter uncertainty on the prediction of two more complex flames: a 2D premixed flame kernel and the ignition of a hydrogen jet issuing into a heated chamber. The former represents a combustion regime similar to the target experiments used to calibrate the mechanism and the latter represents a different combustion regime. For the premixed flame, the net amount of product after a given time interval has a standard deviation of less than 2% whereas the standard deviation of the ignition time for the jet is more than 10%. The samples used for these studies are posted online. These results indicate the degree to which parameters consistent with the target experiments constrain predicted behavior in different combustion regimes. This process provides a framework for both identifying reactions for further study from candidate mechanisms as well as combining uncertainty quantification and propagation to, ultimately, tie uncertainty in laboratory flame experiments to uncertainty in end-use numerical predictions of more complicated scenarios.

In this paper, we explore the suitability of upcoming novel computingtechnologies, in particular adiabatic annealing based quantum computers, tosolve fluid dynamics problems that form a critical component of several scienceand engineering applications. We start with simple flows with well-studied flowproperties, and provide a framework to convert such systems to a form amenablefor deployment on such quantum annealers. We analyze the solutions obtainedboth qualitatively and quantitatively as well as the sensitivities of thevarious solution selection schemes on the obtained solution.

This work introduces a general framework for constructing high-order, linearly stable, partitioned solvers for multiphysics problems from a monolithic implicit–explicit Runge–Kutta (IMEX-RK) discretization of the semi-discrete equations. The generic multiphysics problem is modeled as a system of n systems of partial differential equations where the ith subsystem is coupled to the other subsystems through a coupling term that can depend on the state of all the other subsystems. This coupled system of partial differential equations reduces to a coupled system of ordinary differential equations via the method of lines where an appropriate spatial discretization is applied to each subsystem. The coupled system of ordinary differential equations is taken as a monolithic system and discretized using an IMEX-RK discretization with a specific implicit–explicit decomposition that introduces the concept of a predictor for the coupling term. We propose four coupling predictors that enable the monolithic system to be solved in a partitioned manner, i.e., subsystem-by-subsystem, and preserve the IMEX-RK structure and therefore the design order of accuracy of the monolithic scheme. The four partitioned solvers that result from these predictors are high-order accurate, allow for maximum re-use of existing single-physics software, and two of the four solvers allow the subsystems to be solved in parallel at a given stage and time step. We also analyze the stability of a coupled, linear model problem with a specific coupling structure and show that one of the partitioned solvers achieves unconditional linear stability for this problem, while the others are unconditionally stable only for certain values of the coupling strength. We demonstrate the performance of the proposed partitioned solvers on several classes of multiphysics problems including a simple linear system of ODEs, advection–diffusion–reaction systems, fluid–structure interaction problems, and particle-laden flows, where we verify the design order of the IMEX schemes and study various stability properties.

We investigate how ideas from covariance localization in numerical weather prediction can be used in Markov chain Monte Carlo (MCMC) sampling of high-dimensional posterior distributions arising in Bayesian inverse problems. To localize an inverse problem is to enforce an anticipated “local” structure by (i) neglecting small off-diagonal elements of the prior precision and covariance matrices; and (ii) restricting the influence of observations to their neighborhood. For linear problems we can specify the conditions under which posterior moments of the localized problem are close to those of the original problem. We explain physical interpretations of our assumptions about local structure and discuss the notion of high dimensionality in local problems, which is different from the usual notion of high dimensionality in function space MCMC. The Gibbs sampler is a natural choice of MCMC algorithm for localized inverse problems and we demonstrate that its convergence rate is independent of dimension for localized linear problems. Nonlinear problems can also be tackled efficiently by localization and, as a simple illustration of these ideas, we present a localized Metropolis-within-Gibbs sampler. Several linear and nonlinear numerical examples illustrate localization in the context of MCMC samplers for inverse problems.

We propose a new method for computing Dynamic Mode Decomposition (DMD)evolution matrices, which we use to analyze dynamical systems. Unlike themajority of existing methods, our approach is based on a variationalformulation consisting of data alignment penalty terms and constitutiveorthogonality constraints. Our method does not make any assumptions on thestructure of the data or their size, and thus it is applicable to a wide rangeof problems including non-linear scenarios or extremely small observation sets.In addition, our technique is robust to noise that is independent of thedynamics and it does not require input data to be sequential. Our key idea isto introduce a regularization term for the forward and backward dynamics. Theobtained minimization problem is solved efficiently using the AlternatingMethod of Multipliers (ADMM) which requires two Sylvester equation solves periteration. Our numerical scheme converges empirically and is similar to aprovably convergent ADMM scheme. We compare our approach to variousstate-of-the-art methods on several benchmark dynamical systems.

An efficient hp-multigrid scheme is presented for local discontinuous Galerkin (LDG) discretizations of elliptic problems, formulated around the idea of separately coarsening the underlying discrete gradient and divergence operators. We show that traditional multigrid coarsening of the primal formulation leads to poor and suboptimal multigrid performance, whereas coarsening of the flux formulation leads to essentially optimal convergence and is equivalent to a purely geometric multigrid method. The resulting operator-coarsening schemes do not require the entire mesh hierarchy to be explicitly built, thereby obviating the need to compute quadrature rules, lifting operators, and other mesh-related quantities on coarse meshes. We show that good multigrid convergence rates are achieved in a variety of numerical tests on two-dimensional and three-dimensional uniform and adaptive Cartesian grids, as well as for curved domains using implicitly defined meshes and for multiphase elliptic interface problems with complex geometry. Extensions, e.g., to non-LDG discretizations and fully unstructured meshes, are briefly discussed.