Publications

Well-conditioned boundary integral methods for the solution of elliptic boundary value problems (BVPs) are powerful tools for static and dynamic physical simulations. When there are many close-to-touching boundaries (e.g., in complex fluids) or when the solution is needed in the bulk, nearly singular integrals must be evaluated at many targets. We show that precomputing a linear map from surface density to an effective source representation renders this task highly efficient, in the common case where each object is “simple”, i.e., its smooth boundary needs only moderately many nodes. We present a kernel-independent method needing only an upsampled smooth surface quadrature, and one dense factorization, for each distinct shape. No (near-)singular quadrature rules are needed. The resulting effective sources are drop-in compatible with fast algorithms, with no local corrections nor bookkeeping. Our extensive numerical tests include 2D FMM-based Helmholtz and Stokes BVPs with up to 1000 objects (281000 unknowns), and a 3D Laplace BVP with 10 ellipsoids separated by 1/30 of a diameter. We include a rigorous analysis for analytic data in 2D and 3D.

Single particle cryo-electron microscopy has become a critical tool in structural biology over the last decade, able to achieve atomic scale resolution in three dimensional models from hundreds of thousands of (noisy) two-dimensional projection views of particles frozen at unknown orientations. This is accomplished by using a suite of software tools to (i) identify particles in large micrographs, (ii) obtain low-resolution reconstructions, (iii) refine those low-resolution structures, and (iv) finally match the obtained electron scattering density to the constituent atoms that make up the macromolecule or macromolecular complex of interest. Here, we focus on the second stage of the reconstruction pipeline: obtaining a low resolution model from picked particle images. Our goal is to create an algorithm that is capable of ab initio reconstruction from small data sets (on the order of a few thousand selected particles). More precisely, we seek an algorithm that is robust, automatic, able to assess particle quality, and fast enough that it can potentially be used to assist in the assessment of the data being generated while the microscopy experiment is still underway.

State redistribution is an algorithm that stabilizes cut cells for embedded boundary grid methods. This work extends the earlier algorithm in several important ways. First, state redistribution is extended to three spatial dimensions. Second, we discuss several algorithmic changes and improvements motivated by the more complicated cut cell geometries that can occur in higher dimensions. In particular, we introduce a weighted version with less dissipation in an easily generalizable framework. Third, we demonstrate that state redistribution can also stabilize a solution update that includes both advective and diffusive contributions. The stabilization algorithm is shown to be effective for incompressible as well as compressible reacting flows. Finally, we discuss the implementation of the algorithm for several exascale-ready simulation codes based on AMReX, demonstrating ease of use in combination with domain decomposition, hybrid parallelism and complex physics.

Simulations with adaptive time-dependent bias enable an efficient exploration of the conformational space of a system. However, the dynamic information is altered by the bias. Infrequent metadynamics recovers the transition rate of crossing a barrier, if the collective variables are ideal and there is no bias deposition near the transition state. Unfortunately, these conditions are not always fulfilled. To overcome these limitations, and inspired by single-molecule force spectroscopy, we use Kramers’ theory for calculating the barrier-crossing rate when a time-dependent bias is added to the system. We assess the efficiency of collective variables parameter by measuring how efficiently the bias accelerates the transitions. We present approximate analytical expressions of the survival probability, reproducing the barrier-crossing time statistics and enabling the extraction of the unbiased transition rate even for challenging cases. We explore the limits of our method and provide convergence criteria to assess its validity.

We introduce a numerical method for the solution of the time-dependent Schrödinger equation with a smooth potential, based on its reformulation as a Volterra integral equation. We present versions of the method both for periodic boundary conditions, and for free space problems with compactly supported initial data and potential. A spatially uniform electric field may be included, making the solver applicable to simulations of light-matter interaction. The primary computational challenge in using the Volterra formulation is the application of a spacetime history dependent integral operator. This may be accomplished by projecting the solution onto a set of Fourier modes, and updating their coefficients from one time step to the next by a simple recurrence. In the periodic case, the modes are those of the usual Fourier series, and the fast Fourier transform (FFT) is used to alternate between physical and frequency domain grids. In the free space case, the oscillatory behavior of the spectral Green's function leads us to use a set of complex-frequency Fourier modes obtained by discretizing a contour deformation of the inverse Fourier transform, and we develop a corresponding fast transform based on the FFT. Our approach is related to pseudospectral methods, but applied to an integral rather than the usual differential formulation. This has several advantages: it avoids the need for artificial boundary conditions, admits simple, inexpensive, high-order implicit time marching schemes, and naturally includes time-dependent potentials. We present examples in one and two dimensions showing spectral accuracy in space and eighth-order accuracy in time for both periodic and free space problems.

Partial differential equations (PDEs) play a dominant role in the mathematical modeling of many complex dynamical processes. Solving these PDEs often requires prohibitively high computational costs, especially when multiple evaluations must be made for different parameters or conditions. After training, neural operators can provide PDEs solutions significantly faster than traditional PDE solvers. In this work, invariance properties and computational complexity of two neural operators are examined for transport PDE of a scalar quantity. Neural operator based on graph kernel network (GKN) operates on graph-structured data to incorporate nonlocal dependencies. Here we propose a modified formulation of GKN to achieve frame invariance. Vector cloud neural network (VCNN) is an alternate neural operator with embedded frame invariance which operates on point cloud data. GKN-based neural operator demonstrates slightly better predictive performance compared to VCNN. However, GKN requires an excessively high computational cost that increases quadratically with the increasing number of discretized objects as compared to a linear increase for VCNN.

The long-term goal of this work is the development of high-fidelity simulation tools for dispersive tsunami propagation. A dispersive model is especially important for short wavelength phenomena such as an asteroid impact into the ocean, and is also important in modeling other events where the simpler shallow water equations are insufficient. Adaptive simulations are crucial to bridge the scales from deep ocean to inundation, but have difficulties with the implicit system of equations that results from dispersive models. We propose a fractional step scheme that advances the solution on separate patches with different spatial resolutions and time steps. We show a simulation with 7 levels of adaptive meshes and onshore inundation resulting from a simulated asteroid impact off the coast of Washington. Finally, we discuss a number of open research questions that need to be resolved for high quality simulations.

With the funding provided by this award, we developed numerical codes for the study of magnetically confined plasmas for fusion applications. Accordingly, our work can be divided into two separate categories: 1) the design and analysis of novel numerical methods providing high accuracy and high efficiency; 2) the study of the equilibrium and stability of magnetically confined plasmas with some of these numerical codes, as well as the study of the nature of the turbulent behavior which may arise in the presence of instabilities. We first developed new numerical schemes based on integral equation methods for the computation of steady-state magnetic configurations in fusion experiments, providing high accuracy for the magnetic field and its derivatives, which are required for stability and turbulence calculations. We employed different integral formulations depending on the application of interest: axisymmetric or non-axisymmetric equilibria, force-free or magnetohydrodynamic equilibria, fixed-boundary equilibria or free-boundary equilibria. While efficient, these methods do not yet apply to plasma boundaries which are not smooth, a situation which is fairly common in magnetic confinement experiments. To address this temporary weakness, we also constructed a new steady-state solver based on the Hybridizable Discontinuous Galerkin (HDG) method, which provides full geometric flexibility. In addition to these numerical tools focused on steady-states, we also contributed to the improvement of the speed and accuracy of codes simulating the plasma dynamics of fusion plasmas, by developing a novel velocity space representation for the efficient solution of kinetic equations, which most accurately describe the time evolution of hot plasmas in fusion experiments. Using the tools discussed above, we studied several questions pertaining to the equilibrium and stability of magnetically confined plasmas. In particular, we derived a new simple model for axisymmetric devices called tokamaks, to predict how elongated a fusion plasma can be before it becomes unstable and collapses. We also looked at the effect of the shape of the outer plasma surface on key properties of the steady-state magnetic configurations, and how these properties impact turbulence in fusion plasmas, and the corresponding transport of momentum. Likewise, we studied the role of large localized flows on the steady-state magnetic configurations, and how they may influence plasma stability and turbulence. Non-axisymmetric steady-state magnetic configurations are inherently more complex than axisymmetric steady-state configurations, and the subject of ongoing controversies regarding the regularity of the equations determining such steady-states, and their solutions. Implementing an existing NYU code in a new geometry, we studied the nature of the singularity of the solutions observed in the code, and methods to eliminate them. Our main conclusion is that by appropriately tailoring the plasma boundary, it is possible to eliminate the singularities otherwise appearing in our simulations, and to obtain steady-states which appear to be smooth. To gain further insights on incompletely understood turbulence phenomena, we proposed a new reduced model capturing most of these phenomena, which is simple enough to not require expensive numerical simulations on massive supercomputers to investigate them. We demonstrated the strong similarity between our simulations and published results obtained from computationally expensive simulations, and plan to rely on our reduced model to identify the key mechanisms determining the evolution and strength turbulent driven transport in fusion plasmas. Finally, we proposed a new framework for tokamak reactor design studies, enabling us to consider the relative merits of steady-state versus pulsed fusion reactors. We found that pulsed fusion reactors may benefit most from recent advances in magnet technology, and the availability of very high field magnets. As such, they may become more desirable than steady-state tokamak reactors for cost efficient electricity generation.

A new scheme is proposed to construct a $\mathcal{C}^n\$ function extension for smooth functions defined on a smooth domain $D\in \mathbb{R}^d$. Unlike the PUX scheme, which requires the extrapolation of the volume grid via an expensive ill-conditioned least squares fitting, the scheme relies on an explicit formula consisting of a linear combination of function values in $D,$ which only extends the function along the normal direction. To be more precise, the $\mathcal{C}^n$ extension requires only $n+1$ function values along the normal directions in the original domain and ensures $\mathcal{C}^n$ smoothness by construction. When combined with a shrinking function and a smooth window function, the scheme can be made stable and robust for a broad class of domains with complex smooth boundary.

In this paper, the random iterative method is introduced to massive multiple-input multiple-output (MIMO) systems for the efficient downlink linear precoding. By adopting the random sampling into the traditional iterative methods, the matrix inversion within the linear precoding schemes can be approximated statistically, which not only achieves a faster exponential convergence with low complexity but also experiences a global convergence without suffering from the various convergence requirements. Specifically, based on the random iterative method, the randomized iterative precoding algorithm (RIPA) is firstly proposed and we show its approximation error decays exponentially and globally along with the number of iterations. Then, with respect to the derived convergence rate, the concept of conditional sampling is introduced, so that further optimization and enhancement are carried out to improve both the convergence and the efficiency of the randomized iterations. After that, based on the equivalent iteration transformation, the modified randomized iterative precoding algorithm (MRIPA) is presented, which achieves a better precoding performance with low-complexity for various scenarios of massive MIMO. Finally, simulation results based on downlink precoding in massive MIMO systems are given to show the system gains of RIPA and MRIPA in terms of performance and complexity.