Manas Rachh

Eduardo Corona, Leslie Greengard, Manas Rachh et al.

We present a new derivation of a boundary integral equation (BIE) for simulating the three-dimensional dynamics of arbitrarily-shaped rigid particles of genus zero immersed in a Stokes fluid, on which are prescribed forces and torques. Our method is based on a single-layer representation and leads to a simple second-kind integral equation. It avoids the use of auxiliary sources within each particle that play a role in some classical formulations. We use a spectrally accurate quadrature scheme to evaluate the corresponding layer potentials, so that only a small number of spatial discretization points per particle are required. The resulting discrete sums are computed in $\mathcal{O}(n)$ time, where $n$ denotes the number of particles, using the fast multipole method (FMM). The particle positions and orientations are updated by a high-order time-stepping scheme. We illustrate the accuracy, conditioning and scaling of our solvers with several numerical examples.

Yuanxun Bao, Manas Rachh, Eric Keaveny et al.

We present a fluctuating boundary integral method (FBIM) for overdamped Brownian Dynamics (BD) of two-dimensional periodic suspensions of rigid particles of complex shape immersed in a Stokes fluid. We develop a novel approach for generating Brownian displacements that arise in response to the thermal fluctuations in the fluid. Our approach relies on a first-kind boundary integral formulation of a mobility problem in which a random surface velocity is prescribed on the particle surface, with zero mean and covariance proportional to the Green's function for Stokes flow (Stokeslet). This approach yields an algorithm that scales linearly in the number of particles for both deterministic and stochastic dynamics, handles particles of complex shape, achieves high order of accuracy, and can be generalized to three dimensions and other boundary conditions. We show that Brownian displacements generated by our method obey the discrete fluctuation-dissipation balance relation (DFDB). Based on a recently-developed Positively Split Ewald method [A. M. Fiore, F. Balboa Usabiaga, A. Donev and J. W. Swan, J. Chem. Phys., 146, 124116, 2017], near-field contributions to the Brownian displacements are efficiently approximated by iterative methods in real space, while far-field contributions are rapidly generated by fast Fourier-space methods based on fluctuating hydrodynamics. FBIM provides the key ingredient for time integration of the overdamped Langevin equations for Brownian suspensions of rigid particles. We demonstrate that FBIM obeys DFDB by performing equilibrium BD simulations of suspensions of starfish-shaped bodies using a random finite difference temporal integrator.

Mo Zhou, Jiequn Han, Manas Rachh et al.

We consider the inverse acoustic obstacle problem for sound-soft star-shaped obstacles in two dimensions wherein the boundary of the obstacle is determined from measurements of the scattered field at a collection of receivers outside the object. One of the standard approaches for solving this problem is to reformulate it as an optimization problem: finding the boundary of the domain that minimizes the $L^2$ distance between computed values of the scattered field and the given measurement data. The optimization problem is computationally challenging since the local set of convexity shrinks with increasing frequency and results in an increasing number of local minima in the vicinity of the true solution. In many practical experimental settings, low frequency measurements are unavailable due to limitations of the experimental setup or the sensors used for measurement. Thus, obtaining a good initial guess for the optimization problem plays a vital role in this environment. We present a neural network warm-start approach for solving the inverse scattering problem, where an initial guess for the optimization problem is obtained using a trained neural network. We demonstrate the effectiveness of our method with several numerical examples. For high frequency problems, this approach outperforms traditional iterative methods such as Gauss-Newton initialized without any prior (i.e., initialized using a unit circle), or initialized using the solution of a direct method such as the linear sampling method. The algorithm remains robust to noise in the scattered field measurements and also converges to the true solution for limited aperture data. However, the number of training samples required to train the neural network scales exponentially in frequency and the complexity of the obstacles considered. We conclude with a discussion of this phenomenon and potential directions for future research.

Manas Rachh, Andreas Klöckner, Michael O'Neil

The use of integral equation methods for the efficient numerical solution of PDE boundary value problems requires two main tools: quadrature rules for the evaluation of layer potential integral operators with singular kernels, and fast algorithms for solving the resulting dense linear systems. Classically, these tools were developed separately. In this work, we present a unified numerical scheme based on coupling Quadrature by Expansion, a recent quadrature method, to a customized Fast Multipole Method (FMM) for the Helmholtz equation in two dimensions. The method allows the evaluation of layer potentials in linear-time complexity, anywhere in space, with a uniform, user-chosen level of accuracy as a black-box computational method. Providing this capability requires geometric and algorithmic considerations beyond the needs of standard FMMs as well as careful consideration of the accuracy of multipole translations. We illustrate the speed and accuracy of our method with various numerical examples. Keywords: Layer Potentials; Singular Integrals; Quadrature; High-order accuracy; Integral equations; Helmholtz equation; Fast multipole method.

Florian Pausinger, Manas Rachh, Stefan Steinerberger

Jittered Sampling is a refinement of the classical Monte Carlo sampling method. Instead of picking $n$ points randomly from $[0,1]^2$, one partitions the unit square into $n$ regions of equal measure and then chooses a point randomly from each partition. Currently, no good rules for how to partition the space are available. In this paper, we present a solution for the special case of subdividing the unit square by a decreasing function into two regions so as to minimize the expected squared $\mathcal{L}_2-$discrepancy. The optimal partitions are given by a \textit{highly} nonlinear integral equation for which we determine an approximate solution. In particular, there is a break of symmetry and the optimal partition is not into two sets of equal measure. We hope this stimulates further interest in the construction of good partitions.

Jeremy G. Hoskins, Alan E. Lindsay, Manas Rachh

In this paper we present a new boundary integral equation formulation for the solution of the elastostatic traction boundary value problem in two and three dimensions. The approach relies on the introduction of new layer potentials, called string kernels, which are based on modifications of the Boussinesq-Cerruti family of half-space solutions. We prove that the resulting integral equations are second-kind integral equations, and show that they are well-behaved in the incompressible limit. We illustrate the performance of the method with several numerical examples.

Jeffrey Galkowski, Manas Rachh, Euan A. Spence

This paper is concerned with solving the Helmholtz exterior Dirichlet and Neumann problems with large wavenumber $k$ and smooth obstacles using the standard second-kind boundary integral equations. We consider Galerkin and collocation methods -- with subspaces consisting of $\textit{either}$ piecewise polynomials (in 2-d for collocation, in any dimension for Galerkin) $\textit{or}$ trigonometric polynomials (in 2-d) -- as well as a fully discrete quadrature (Nyström) method based on trigonometric polynomials (in 2-d). For each of these methods, we prove -- in many cases for the first time -- rigorous results about the fundamental question: how quickly must the number of degrees of freedom (the dimension of the approximation space) grow with $k$ to maintain accuracy of the computed solution? Importantly, we determine which of these methods suffer from $\textit{the pollution effect}$. That is, we address the question: must the number of points per wavelength $\to \infty$ to maintain accuracy as $k\to\infty$?

Guillaume Bal, Jeremy Hoskins, Solomon Quinn et al.

This paper concerns a boundary integral formulation for the two-dimensional massive Dirac equation. The mass term is assumed to jump across a one-dimensional interface, which models a transition between two insulating materials. This jump induces surface waves that propagate outward along the interface but decay exponentially in the transverse direction. After providing a derivation of our integral equation, we prove that it has a unique solution for almost all choices of parameters using holomorphic perturbation theory. We then extend these results to a Dirac equation with two interfaces. Finally, we implement a fast numerical method for solving our boundary integral equations and present several numerical examples of solutions and scattering effects.

Zijian Wang, Jeremy G. Hoskins, Manas Rachh et al.

In this paper, we study the scale-invariant quantity \[\mathcal{G}(Ω)=\frac{\|\partial_n u_1\|_{L^\infty(\partialΩ)}}{λ_1},\]where $u_1$ is the first $L^2$-normalized Dirichlet Laplace eigenfunction of a Euclidean domain $Ω$ and $λ_1$ is its eigenvalue. This is related to the peak boundary heat flux in the long time limit. For convex domains we prove that $\|\partial_n u_1\|_{L^\infty(\partialΩ)}$ is upper-bounded by a (domain-independent) constant multiple of $λ_1$. Using layer potentials, we derive shape-derivative formulae for efficient gradient computations. When combined with high-order Nyström discretization, a fast boundary integral equation solver, and eigenvalue rootfinding, this allows us to numerically optimize $\mathcal{G}$ over a class of rounded polygonal discretized domains. Based on extensive numerical experiments, we then conjecture that, over the set of convex domains, $\mathcal{G}$ is maximized by the semidisk, with the peak flux at the center of the diameter. To lend analytical support to this conjecture, we prove that the semidisk is a critical point of $\mathcal{G}$ under infinitesimal perturbations of its circular arc.

Travis Askham, Manas Rachh

The eigenvalues and eigenfunctions of the Stokes operator have been the subject of intense analytical investigation and have applications in the study and simulation of the Navier-Stokes equations. As the Stokes operator is a fourth-order operator, computing these eigenvalues and the corresponding eigenfunctions is a challenging task, particularly in complex geometries and at high frequencies. The boundary integral equation (BIE) framework provides robust and scalable eigenvalue computations due to (a) the reduction in the dimension of the problem to be discretized and (b) the absence of high frequency "pollution" when using a Green's function to represent propagating waves. In this paper, we detail the theoretical justification for a BIE approach to the Stokes eigenvalue problem on simply and multiply-connected planar domains, which entails a treatment of the uniqueness theory for oscillatory Stokes equations on exterior domains. Then, using well-established techniques for discretizing BIEs, we present numerical results which confirm the analytical claims of the paper and demonstrate the efficiency of the overall approach.

George C. Linderman, Manas Rachh, Jeremy G. Hoskins et al.

t-distributed Stochastic Neighborhood Embedding (t-SNE) is a method for dimensionality reduction and visualization that has become widely popular in recent years. Efficient implementations of t-SNE are available, but they scale poorly to datasets with hundreds of thousands to millions of high dimensional data-points. We present Fast Fourier Transform-accelerated Interpolation-based t-SNE (FIt-SNE), which dramatically accelerates the computation of t-SNE. The most time-consuming step of t-SNE is a convolution that we accelerate by interpolating onto an equispaced grid and subsequently using the fast Fourier transform to perform the convolution. We also optimize the computation of input similarities in high dimensions using multi-threaded approximate nearest neighbors. We further present a modification to t-SNE called "late exaggeration," which allows for easier identification of clusters in t-SNE embeddings. Finally, for datasets that cannot be loaded into the memory, we present out-of-core randomized principal component analysis (oocPCA), so that the top principal components of a dataset can be computed without ever fully loading the matrix, hence allowing for t-SNE of large datasets to be computed on resource-limited machines.

Manas Rachh, Travis Askham

We present a novel integral representation for the biharmonic Dirichlet problem. To obtain the representation, the Dirichlet problem is first converted into a related Stokes problem for which the Sherman-Lauricella integral representation can be used. Not all potentials for the Dirichlet problem correspond to a potential for Stokes flow, and vice-versa, but we show that the integral representation can be augmented and modified to handle either simply or multiply connected domains. The resulting integral representation has a kernel which behaves better on domains with high curvature than existing representations. Thus, this representation results in more robust computational methods for the solution of the Dirichlet problem of the biharmonic equation and we demonstrate this with several numerical examples.

Xiuyuan Cheng, Manas Rachh, Stefan Steinerberger

We study directed, weighted graphs $G=(V,E)$ and consider the (not necessarily symmetric) averaging operator $$ (\mathcal{L}u)(i) = -\sum_{j \sim_{} i}{p_{ij} (u(j) - u(i))},$$ where $p_{ij}$ are normalized edge weights. Given a vertex $i \in V$, we define the diffusion distance to a set $B \subset V$ as the smallest number of steps $d_{B}(i) \in \mathbb{N}$ required for half of all random walks started in $i$ and moving randomly with respect to the weights $p_{ij}$ to visit $B$ within $d_{B}(i)$ steps. Our main result is that the eigenfunctions interact nicely with this notion of distance. In particular, if $u$ satisfies $\mathcal{L}u = λu$ on $V$ and $$ B = \left\{ i \in V: - \varepsilon \leq u(i) \leq \varepsilon \right\} \neq \emptyset,$$ then, for all $i \in V$, $$ d_{B}(i) \log{\left( \frac{1}{|1-λ|} \right) } \geq \log{\left( \frac{ |u(i)| }{\|u\|_{L^{\infty}}} \right)} - \log{\left(\frac{1}{2} + \varepsilon\right)}.$$ $d_B(i)$ is a remarkably good approximation of $|u|$ in the sense of having very high correlation. The result implies that the classical one-dimensional spectral embedding preserves particular aspects of geometry in the presence of clustered data. We also give a continuous variant of the result which has a connection to the hot spots conjecture.