High-order close evaluation of Laplace layer potentials: A differential geometric approach
H. Zhu and S. Veerapaneni
Submitted, 2021.

This paper presents a new approach for solving the close evaluation problem in three dimensions, commonly encountered while solving linear elliptic partial differential equations via potential theory. The goal is to evaluate layer potentials close to the boundary over which they are defined. The approach introduced here converts these nearly-singular integrals on a patch of the boundary to a set of non-singular line integrals on the patch boundary using the Stokes theorem on manifolds. A function approximation scheme based on harmonic polynomials is designed to express the integrand in a form that is suitable for applying the Stokes theorem. As long as the data—the boundary and the density function—is given in a high-order format, the double-layer potential and its derivatives can be evaluated with high-order accuracy using this scheme both on and off the boundary. In particular, we present numerical results demonstrating seventh-order convergence on a smooth, warped torus example achieving 10-digit accuracy in evaluating double layer potential at targets that are arbitrarily close to the boundary.

Fast and accurate solvers for simulating Janus particle suspensions in Stokes flow
R. Kohl, E. Corona, V. Cheruvu and S. Veerapaneni
Submitted, 2021.

We present a novel computational framework for simulating suspensions of rigid spherical Janus particles in Stokes flow. We show that long-range Janus particle interactions for a wide array of applications may be resolved using fast, spectrally accurate boundary integral methods tailored to polydisperse suspensions of spherical particles. These are incorporated into our rigid body Stokes platform. Our approach features the use of spherical harmonic expansions for spectrally accurate integral operator evaluation, complementarity-based collision resolution, and optimal O(n) scaling with the number of particles when accelerated via fast summation techniques. We demonstrate the flexibility of our platform through three key examples of Janus particle systems prominent in biomedical applications: amphiphilic, bipolar electric and phoretic particles. We formulate Janus particle interactions in boundary integral form and showcase characteristic self-assembly and complex collective behavior for each particle type.

Optimal ciliary locomotion of axisymmetric microswimmers
H. Guo, H. Zhu, R. Liu, M. Bonnet and S. Veerapaneni
Journal of Fluid Mechanics, 2021.

Many biological microswimmers locomote by periodically beating the densely-packed cilia on their cell surface in a wave-like fashion. While the swimming mechanisms of ciliated microswimmers have been extensively studied both from the analytical and the numerical point of view, the optimization of the ciliary motion of microswimmers has received limited attention, especially for non-spherical shapes. In this paper, using an envelope model for the microswimmer, we numerically optimize the ciliary motion of a ciliate with an arbitrary axisymmetric shape. The forward solutions are found using a fast boundary integral method, and the efficiency sensitivities are derived using an adjoint-based method. Our results show that a prolate microswimmer with a 2:1 aspect ratio shares similar optimal ciliary motion as the spherical microswimmer, yet the swimming efficiency can increase two-fold.

Optimal slip velocities of micro-swimmers with arbitrary axisymmetric shapes
H. Guo, H. Zhu, R. Liu, M. Bonnet and S. Veerapaneni
Journal of Fluid Mechanics, Vol. 910, 2021.

This article presents a computational approach for determining the optimal slip velocities on any given shape of an axisymmetric micro-swimmer suspended in a viscous fluid. The objective is to minimize the power loss to maintain a target swimming speed, or equivalently to maximize the efficiency of the micro-swimmer. Owing to the linearity of the Stokes equations governing the fluid motion, we show that this PDE-constrained optimization problem reduces to a simpler quadratic optimization problem, whose solution is found using a high-order accurate boundary integral method. We consider various families of shapes parameterized by the reduced volume and compute their swimming efficiency. Among those, prolate spheroids were found to be the most efficient micro-swimmer shapes for a given reduced volume. We propose a simple shape-based scalar metric that can determine whether the optimal slip on a given shape makes it a pusher, a puller or a neutral swimmer.

Natural evolution strategies and variational Monte Carlo
T. Zhao, G. Carleo, J.Stokes and S. Veerapaneni
Machine Learning: Science and Technology, Vol. 2, 2021.

A notion of quantum natural evolution strategies is introduced, which provides a geometric synthesis of a number of known quantum/classical algorithms for performing classical black-box optimization. Recent work of Gomes et al. [2019] on combinatorial optimization using neural quantum states is pedagogically reviewed in this context, emphasizing the connection with natural evolution strategies. The algorithmic framework is illustrated for approximate combinatorial optimization problems, and a systematic strategy is found for improving the approximation ratios. In particular it is found that natural evolution strategies can achieve state-of-art approximation ratios for Max-Cut, at the expense of increased computation time.

An active learning framework for constructing high-fidelity mobility maps
G. Marple, D. Gorsich, P. Jayakumar and S. Veerapaneni
IEEE Transactions on Vehicular Technology, 2021.

Recent work at the U.S. Army CCDC Ground Vehicle Systems Center has shown that machine learning classifiers can quickly construct high-fidelity mobility maps. Training these classifiers, on the other hand, is still a challenge, since each data instance is labeled by performing a computationally intensive, physics-based simulation. In this paper we introduce an active learning framework, based on the query-by-bagging algorithm, that substantially reduces the number of simulations needed to train a classifier. Experimental results suggest that our sampling algorithm can train a neural network, with higher accuracy, using less than half the number of simulations when compared to random sampling.

Simulating cilia-driven mixing and transport in complex geometries
H. Guo, H. Zhu and S. Veerapaneni
Physical Review Fluids, Vol. 5, 2020.

Cilia and flagella are self-actuated microtubule-based structures that are present on many cell surfaces, ranging from the outer surface of single-cell organisms to the internal epithelial surfaces in larger animals. A fast and robust numerical method that simulates the coupled hydrodynamics of cilia and the constituent particles in the fluid such as rigid particles, drops or cells would be useful to not only understand several disease and developmental pathologies due to ciliary dysfunction but also to design microfluidic chips with ciliated cultures for some targeted functionality -- e.g., maximizing fluid transport or particle mixing. In this paper, we develop a hybrid numerical method that employs a boundary integral method for handling the confining geometries and the constituent rigid particles and the method of regularized Stokeslets for handling the cilia. We provide several examples demonstrating the effects of confining geometries on cilia-generated fluid mixing as well as the cilia-particle hydrodynamics.

Fast integral equation methods for linear and semilinear heat equations in moving domains
J. Wang, L. Greengard, S. Jiang and S. Veerapaneni
Submitted, 2020.

We present a family of integral equation-based solvers for the linear or semilinear heat equation in complicated moving (or stationary) geometries. This approach has significant advantages over more standard finite element or finite difference methods in terms of accuracy, stability and space-time adaptivity. In order to be practical, however, a number of technical capabilites are required: fast algorithms for the evaluation of heat potentials, high-order accurate quadratures for singular and weakly integrals over space-time domains, and robust automatic mesh refinement and coarsening capabilities. We describe all of these components and illustrate the performance of the approach with numerical examples in two space dimensions.

A scalable computational platform for particulate Stokes suspensions
W. Yan, E. Corona, D. Malhotra, S. Veerapaneni and M. Shelley
Journal of Computational Physics, Vol. 416, 2020.

We describe a computational framework for simulating suspensions of rigid particles in Newtonian Stokes flow. One central building block is a collision-resolution algorithm that overcomes the numerical constraints arising from particle collisions. This algorithm extends the well-known complementarity method for non-smooth multi-body dynamics to resolve collisions in dense rigid body suspensions. This approach formulates the collision resolution problem as a linear complementarity problem with geometric ‘non-overlapping’ constraints imposed at each timestep. It is then reformulated as a constrained quadratic programming problem and the Barzilai-Borwein projected gradient descent method is applied for its solution... Further, we describe a fast, parallel, and spectrally-accurate boundary integral method tailored for spherical particles, capable of resolving lubrication effects. We show weak and strong parallel scalings up to 8 × 10^4 particles with approximately 4 × 10^7 degrees of freedom on 1792 cores. We demonstrate the versatility of this framework with several examples, including sedimentation of particle clusters, and active matter systems composed of ensembles of particles driven to rotate.

Solution of Stokes flow in complex nonsmooth 2D geometries via a linear-scaling high-order adaptive integral equation scheme
B. Wu, H. Zhu, A. Barnett and S. Veerapaneni
Journal of Computational Physics, Vol. 410, 2020.

We present a fast, high-order accurate and adaptive boundary integral scheme for solving the Stokes equations in complex, possibly nonsmooth, geometries in two dimensions. The key ingredient is a set of panel quadrature rules capable of evaluating weakly-singular, nearly-singular and hyper-singular integrals to high accuracy. Near-singular integral evaluation, in particular, is done using an extension of the scheme developed in J. Helsing and R. Ojala, JCP (2008). The boundary of the given geometry is ``panelized'' automatically to achieve user-prescribed precision. We found that this adaptive mesh refinement procedure works well in practice even in the case of complex geometries with large number of corners. In one example, for instance, a model 2D vascular network with 378 corners required less than 200K discretization points to obtain a 9-digit solution accuracy.

Shape optimization of Stokesian peristaltic pumps using boundary integral methods
M. Bonnet, R. Liu and S. Veerapaneni
Advances in Computational Mathematics, Vol. 46(18), 2020.

This article presents a new boundary integral approach for finding optimal shapes of peristaltic pumps that transport a viscous fluid. Formulas for computing the shape derivatives of the standard cost functionals and constraints are derived. They involve evaluating physical variables (traction, pressure, etc.) on the boundary only. By emplyoing these formulas in conjuction with a boundary integral approach for solving forward and adjoint problems, we completely avoid the issue of volume remeshing when updating the pump shape as the optimization proceeds. This leads to significant cost savings and we demonstrate the performance on several numerical examples.

Scalable Solvers for Cone Complementarity Problems in Frictional Multibody Dynamics
S. De, E. Corona, P. Jayakumar and S. Veerapaneni
Proceedings of IEEE Conference on High Performance Extreme Computing, 2019.

We present an efficient, hybrid MPI/OpenMP framework for the cone complementarity formulation of large-scale rigid body dynamics problems with frictional contact. Data is partitioned among MPI processes using a Morton encoding in order to promote data locality and minimize communication. We parallelize the state-of-the-art first and second-order solvers for the resulting cone complementarity optimization problems. Our approach is highly scalable, enabling the solution of dense, large-scale multibody problems; a sedimentation simulation involving 256 million particles (~324 million contacts on average) was resolved using 512 cores in less than half-hour per time-step.

Hydrodynamics and rheology of a vesicle doublet suspension
B. Quaife, S. Veerapaneni and Y.-N. Young
Physical Review Fluids, Vol. 4(10), 2019. *Editor's suggestion*

The dynamics of an adhesive vesicle doublet under various flow conditions is investigated numerically using a high-order, adaptive-in-time boundary integral method. In a quiescent flow, two nearby vesicles move slowly towards each other under the adhesive potential, pushing out fluid between them to form a vesicle doublet at equilibrium. A lubrication analysis on such draining of a thin film gives the dependencies of draining time on adhesion strength and separation distance that are in good agreement with numerical results. In a microfluid trap where the stagnation of an extensional flow is dynamically placed in the middle of a vesicle doublet through an active control loop, novel dynamics of a vesicle doublet is observed. Numerical simulations show that there exists a critical extensional flow rate above which adhesive interaction is overcome by the diverging stream, thus providing a simple method to measure the adhesion strength between two vesicle membranes...

Tensor train accelerated solvers for nonsmooth rigid body dynamics
E. Corona, D. Gorsich, P. Jayakumar and S. Veerapaneni
Applied Mechanics Reviews, Vol. 71(5), 2019.

...In this work, we propose a novel acceleration for the solution of Newton step linear systems in second order methods using low-rank compression based fast direct solvers, leveraging on recent direct solver techniques for structured linear systems arising from differential and integral equations. We employ the Quantized Tensor Train (QTT) decomposition to produce efficient approximate representations of the system matrix and its inverse...We demonstrate compressibility of the Newton step matrices in Primal Dual Interior Point (PDIP) methods as applied to the multibody dynamics problem. Using a number of numerical tests, we demonstrate that this approach displays sublinear scaling of precomputation costs, may be efficiently updated across Newton iterations as well as across simulation time steps, and leads to a fast, optimal complexity solution of the Newton step.

Electrohydrodynamics of deflated vesicles: budding, rheology and pairwise interactions
B. Wu and S. Veerapaneni
Journal of Fluid Mechanics, Vol. 867, pp. 334-337, 2019.

The electrohydrodynamics of vesicle suspensions is characterized by studying their pairwise interactions in applied DC electric fields in two dimensions. In the dilute limit, the rheology of the suspension is shown to vary nonlinearly with the electric conductivity ratio of the interior and exterior fluids...When two vesicles are initially un-aligned with the external electric field, three different responses are observed when the key parameters are varied: (i) chain formation–they self-assemble to form a chain that is aligned along the field direction, (ii) circulatory motion–they rotate about each other, (iii) oscillatory motion–they form a chain but oscillate about each other.

Boundary integral equation analysis for suspension of spheres in Stokes flow
E. Corona and S. Veerapaneni
Journal of Computational Physics, Vol. 362, pp. 327-345, 2018.

We show that the standard boundary integral operators, defined on the unit sphere, for the Stokes equations diagonalize on a specific set of vector spherical harmonics and provide formulas for their spectra. We also derive analytical expressions for evaluating the operators away from the boundary. When two particle are located close to each other, we use a truncated series expansion to compute the hydrodynamic interaction. On the other hand, we use the standard spectrally accurate quadrature scheme to evaluate smooth integrals on the far-field, and accelerate the resulting discrete sums using the fast multipole method. We employ this discretization scheme to analyze several boundary integral formulations of interest including those arising in porous media flow, active matter and magneto-hydrodynamics of rigid particles. We provide numerical results verifying the accuracy and scaling of their evaluation.

A unified integral equation scheme for doubly-periodic Laplace and Stokes boundary value problems in two dimensions
A. Barnett, G. Marple, S. Veerapaneni and L. Zhao
Communications on Pure and Applied Mathematics, 2018.

We present a spectrally-accurate scheme to turn a boundary integral formulation for an elliptic PDE on a single unit cell geometry into one for the fully periodic problem. Applications include computing the effective permeability of composite media (homogenization), and microfluidic chip design. Our basic idea is to exploit a small least squares solve to apply periodicity without ever handling periodic Green’s functions....The scheme is simple (no lattice sums, Ewald methods, nor particle meshes are required), allows adaptivity, and is essentially dimension- and PDE-independent, so would generalize without fuss to 3D and to other non-oscillatory elliptic problems such as elastostatics. We incorporate recently developed spectral quadratures that accurately handle close-to-touching geometries. We include many numerical examples, and provide a software implementation.

Dynamics of a multicomponent vesicle in shear flow
K. Liu, G. Marple, S. Li, S. Veerapaneni and J. Lowengrub
Soft Matter, Vol. 13, pp. 3521-3531, 2017.

We study the fully nonlinear, nonlocal dynamics of two-dimensional multicomponent vesicles in a shear flow with matched viscosity of the inner and outer fluids. Using a nonstiff, pseudo-spectral boundary integral method, we investigate dynamical patterns induced by inhomogeneous bending for a two phase system. Numerical results reveal that there exist novel phase-treading and tumbling mechanisms that cannot be observed for a homogeneous vesicle. In particular, unlike the well-known steady tank-treading dynamics characterized by a fixed inclination angle, here the phase-treading mechanism leads to unsteady periodic dynamics with an oscillatory inclination angle. When the average phase concentration is around 1/2, we observe tumbling dynamics even for very low shear rate, and the excess length required for tumbling is significantly smaller than the value for the single phase case. We summarize our results in phase diagrams in terms of the excess length, shear rate, and concentration of the soft phase. These findings go beyond the well known dynamical regimes of a homogeneous vesicle and highlight the level of complexity of vesicle dynamics in a fluid due to heterogeneous material properties.

An integral equation formulation for rigid bodies in Stokes flow in three dimensions
E. Corona, L. Greengard, M. Rachh and S. Veerapaneni
Journal of Computational Physics, Vol. 332, pp. 504-519, 2017.

We present a new derivation of a boundary integral equation (BIE) for simulating the threedimensional 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 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.

A fast algorithm for simulating multiphase flows through periodic geometries of arbitrary shape
G. Marple, A. Barnett, A. Gillman and S. Veerapaneni
SIAM Journal on Scientific Computing, Vol. 38(5), pp. B740-B772, 2016.

This paper presents a new boundary integral equation (BIE) method for simulating particulate and multiphase flows through periodic channels of arbitrary smooth shape in two dimensions.... Rather than relying on the periodic Green's function as classical BIE methods do, the method combines the free-space Green's function with a small auxiliary basis, and imposes periodicity as an extra linear condition. As a result, we can exploit existing free-space solver libraries, quadratures, and fast algorithms, and handle a large number of vesicles in a geometrically complex channel.... New constraint-correction formulas are introduced that preserve reduced areas of vesicles, independent of the number of time steps taken.... Numerical experiments illustrate that a simulation with N=128,000 can be evolved in less than a minute per time step on a laptop.

Integral equation methods for vesicle electrohydrodynamics in three dimensions
S. Veerapaneni
Journal of Computational Physics, Vol. 326, pp. 278-289, 2016.

In this paper, we develop a new boundary integral equation formulation that describes the coupled electro- and hydro-dynamics of a vesicle suspended in a viscous fluid and subjected to external flow and electric fields. The dynamics of the vesicle are characterized by a competition between the elastic, electric and viscous forces on its membrane. The classical Taylor-Melcher leaky-dielectric model is employed for the electric response of the vesicle and the Helfrich energy model combined with local inextensibility is employed for its elastic response. The coupled governing equations for the vesicle position and its transmembrane electric potential are solved using a numerical method that is spectrally accurate in space and first-order in time. The method uses a semi-implicit time-stepping scheme to overcome the numerical stiffness associated with the governing equations. Using this method, we simulate, for the first time, the buckling process of vesicles under applied electric field in three dimensions.

Simple and efficient representations for the fundamental solutions of Stokes flow in a half-space
Z. Gimbutas, L. Greengard and S. Veerapaneni
Journal of Fluid Mechanics, Vol. 776, Aug. 2015.

We derive new formulas for the fundamental solutions of slow, viscous flow, governed by the Stokes equations, in a half-space. They are simpler than the classical representations obtained by Blake and collaborators, and can be efficiently implemented using existing fast solvers libraries. We show, for example, that the velocity field induced by a Stokeslet can be annihilated on the boundary (to establish a zero slip condition) using a single reflected Stokeslet combined with a single Papkovich-Neuber potential that involves only a scalar harmonic function. The new representation has a physically intuitive interpretation.

Gating of a mechanosensitive channel due to cellular flows
O. Pak, Y. Young, G. Marple, S. Veerapaneni and H. Stone
Proceedings of the National Academy of Sciences, Vol. 112(32), pp. 9822-9827, 2015.

In this work we construct a multiscale continuum model for a mechanosensitive (MS) channel gated by tension in a lipid bilayer membrane under fluid stress. We illustrate that for typical physiological conditions vesicle hydrodynamics driven by a fluid flow may render a membrane tension sufficiently large to gate a MS channel open. We focus on the dynamic opening/closing of a MS channel in a vesicle membrane under a planar shear flow and a pressure-driven flow across a constriction channel. Our modeling and numerical simulation results quantify the critical flow strength or flow channel geometry for intracellular transport through a MS channel. Furthermore our results suggest that lipid bilayer membrane reconstituted with MS channels can be mechanically gated for permeability under fluid flows that are physiologically relevant and realizable in microfluidic configurations for medical applications. These modeling and simulation results imply that stress-induced intracellular transport across the lipid membrane can be achieved by the gating of reconstituted MS channels, which can be useful for designing drug delivery in medical therapy and understanding complicated mechanotransduction.

Equilibrium shapes of planar elastic membranes
G. Marple, P. Purohit and S. Veerapaneni
Physical Review E, Vol. 92 (1), Jul 2015.

Using a rod theory formulation, we derive equations of state for a thin elastic membrane subjected to two different boundary conditions -- clamped and periodic. The former is applicable to membranes supported on a softer substrate and subjected to uniaxial compression. We show that a wider family of quasi-static equilibrium shapes exist beyond the previously obtained analytical solutions. In the latter case of periodic membranes, we were able to derive exact solutions in terms of elliptic functions. These equilibria are verified by considering a fluid-structure interaction problem of a periodic, length-preserving bilipid membrane modeled by the Helfrich energy immersed in a viscous fluid. Starting from an arbitrary shape, the membrane dynamics to equilibrium are simulated using a boundary integral method.

Boundary integral method for the flow of vesicles with viscosity contrast in three dimensions
A. Rahimian, S. Veerapaneni, D. Zorin and G. Biros
Journal of Computational Physics, Vol. 298, pp. 766-786, 2015.

We propose numerical algorithms for the simulation of the dynamics of three-dimensional vesicles suspended in viscous Stokesian fluid. Our method is an extension of our previous work [S. K. Veerapaneni, A. Rahimian, G. Biros, and D. Zorin. A Fast Algorithm for Simulating Vesicle Flows in Three Dimensions. Journal of Computational Physics, 230(14):5610{5634, 2011] to flows with viscosity contrast. This generalization poses two new types of difficulties: (i) change in the boundary integral formulation of the solution, in which a double-layer Stokes integral is introduced; and (ii) change of the fluid dynamics due to the viscosity contrastof the vesicles.We propose algorithms to address these challenges.

Spectrally-accurate quadratures for evaluation of layer potentials close to the boundary for the 2D Stokes and Laplace equations
A. Barnett, B. Wu and S. Veerapaneni
SIAM Journal on Scientific Computing, Vol. 37(4), pp. B519-B542, 2015.

Dense particulate flow simulations using integral equation methods demand accurate evaluation of Stokes layer potentials on arbitrarily close interfaces. In this paper, we generalize techniques for close evaluationof Laplace double-layer potentials in J. Helsing and R. Ojala, J. Comput. Phys. 227 (2008) 2899–2921. We create a “globally compensated” trapezoid rule quadrature for the Laplace single-layer potential on the interior and exterior of smooth curves. This exploits a complex representation, a product quadrature (in the style of Kress) for the sawtooth function, careful attention to branch cuts, and second-kind barycentric-type formulae for Cauchy integrals and their derivatives. Upon this we build accurate single- and double-layer Stokes potential evaluators by expressing them in terms of Laplace potentials. We test their convergence for vesicle-vesicle interactions, for an extensive set of Laplace and Stokes problems, and when applying the system matrix in a boundary value problem solver in the exterior of multiple close-to-touching ellipses. We achieve typically 12 digits of accuracy using very small numbers of discretization nodes per curve. We provide documented codes for other researchers to use.

Long-wave dynamics of an inextensible planar membrane in an electric field
Y. Young, S. Veerapaneni and M. Miksis
Journal of Fluid Mechanics, Vol. 751, pp. 406-431 , 2014.

We investigate the long-wave nonlinear dynamics of an inextensible capacitive elastic membrane under electric fields. In the lubrication framework we derive a nonlinear equation for the membrane height with an integral constraint. Linear analysis on the tensionless membrane in a dc field gives the linear growth rate in terms of membrane conductance and electric properties in the bulk. The long-wave formulation allows us to analytically derive the equilibrium membrane profile in a dc field. Numerical simulations of an inextensible membrane under ac fields elucidate how variation of the membrane tension correlates to the non-linear membrane dynamics. Different membrane dynamics, such as undulation and flip-flop, is found at different electric field strength and membrane area. In particular a traveling wave on the membrane is found as a response to a periodic ac field in the perpendicular direction.

A fast algorithm for spherical grid rotations and its application to singular quadrature
Z. Gimbutas and S. Veerapaneni
SIAM Journal on Scientific Computing, Vol. 35(6), 2013.

We present a fast and accurate algorithm for evaluating singular integral operators on smooth surfaces that are globally parametrized by
spherical coordinates. Problems of this type arise, for example, in simulating Stokes flows with particulate suspensions and in multi-particle scattering calculations. For smooth surfaces, spherical harmonic expansions are commonly used for geometry representation and the evaluation of the singular integrals is carried out with a spectrally accurate quadrature rule on a set of rotated spherical grids. We propose a new algorithm that interpolates function values on the rotated spherical grids via hybrid nonuniform FFTs. The algorithm is nearly optimal in computational complexity and reduces
the cost of applying the quadrature rule from O(p^5 ) to O(p^4 log p) for a spherical harmonic expansion of degree p.

Integral equation methods for unsteady stokes flow in two dimensions
S. Jiang, S. Veerapaneni and L. Greengard
SIAM Journal on Scientific Computing, Vol. 34(4), 2012.

We present an integral equation formulation for the unsteady Stokes equations in two dimensions. This problem is of interest in its own right, as a model for slow viscous flow, but perhaps more importantly, as an ingredient in the solution of the full, incompressible Navier-Stokes equations. Using the unsteady Green's function, the velocity evolves analytically as a divergence-free vector field. This avoids the need for either the solution of coupled field equations (as in fully implicit PDE-based marching schemes) or the projection of the velocity field onto a divergence free field at each time step (as in operator splitting methods). In addition to discussing the analytic properties of the operators that arise in the integral formulation, we describe a family of high-order accurate numerical schemes and illustrate their performance with several examples.

Dynamics of a compound vesicle in shear flow
S. Veerapaneni, Y.-N. Young, P. M. Vlahovska and J. Blawzdzeiwicz
Physical Review Letters, Vol. 106(15), 2011.

The dynamics of a compound vesicle (a lipid bilayer membrane enclosing a fluid with a suspended particle) in shear flow is investigated using both numerical simulations and theoretical analysis. We find that the non-linear hydrodynamic interaction between the inclusion and the confining membrane gives rise to new features of the vesicle dynamics: Transition from tank--treading to tumbling can occur in the absence of any viscosity mismatch, and a vesicle can swing if the enclosed particle is non-spherical. Our results highlight the complex effects of internal cellular structures have on cell dynamics in microcirculatory flows. For example, parasites in malaria-infected erythrocytes increase cytoplasmic viscosity, which leads to increase in blood viscosity.

A fast algorithm for simulating vesicle flows in three dimensions
S. Veerapaneni, A. Rahimian, G. Biros and D. Zorin
Journal of Computational Physics, Vol. 230(14), 2011.

Vesicles are locally-inextensible fluid membranes that can sustain bending. In this paper, we extend ``A numerical method for simulating the dynamics of 3{D} axisymmetric vesicles suspended in viscous flows'', Veerapaneni et al. Journal of Computational Physics, 228(19), 2009 to general non-axisymmetric vesicle flows in three dimensions.

Although the main components of the algorithm are similar in spirit to the axisymmetric case (spectral approximation in space, semi-implicit time-stepping scheme), important new elements need to be introduced for a full 3D method. In particular, spatial quantities are discretized using spherical harmonics, and quadrature rules for singular surface integrals need to be adapted to this case; an algorithm for surface reparameterization is neeed to ensure sufficient of the time-stepping scheme, and spectral filtering is introduced to maintain reasonable accuracy while minimizing computational costs. To characterize the stability of the scheme and to construct preconditioners for the iterative linear system solvers used in the semi-implicit time-stepping scheme, we perform a spectral analysis of the evolutions operator on the unit sphere.

By introducing these algorithmic components, we obtain a time-stepping scheme that experimentally is unconditionally stable and has a low cost per time step. We present numerical results to analyze the cost and convergence rates of the scheme. To illustrate the applicability of method, we consider a few vesicle-flow interaction problems: a single vesicle in relaxation, sedimentation, and shear flows, and many-vesicle flows.

Petascale direct numerical simulation of blood flow on 200K cores and heterogeneous architectures
A. Rahimian, I. Lashuk, S. Veerapaneni, A. Chandramowlishwaran, D. Malhotra, L. Moon, R. Sampath, A. Shringarpure, J. Vetter, R.Vuduc, D. Zorin and G. Biros
ACM/IEEE Conference on Supercomputing, 2010. *Gordon Bell Prize*

We present a high fidelity numerical simulation of blood flow by directly resolving the interactions of 200 million deformable red blood cells flowing in plasma. This simulation amounts to 90 billion unknowns in space, with numerical experiments typically requiring O(1000) time steps. In terms of the number of cells, we improve the state-of-the art by several orders of magnitude: the previous largest simulation, at the same physical fidelity as ours, resolved the flow of 14 thousand cells. This breakthrough is based on novel algorithms that we designed to enable distributed memory, shared memory, and vectorized/streaming parallelism. We present results on CPU and hybrid CPU-GPU platforms, including the new NVIDIA Fermi architecture and 200,000 cores of ORNL's Jaguar system. For the latter, we achieve over 0.7 Petaflop/s sustained performance. Our work demonstrates the successful simulation of complex phenomena using heterogeneous architectures and programming models at the petascale.

Parallel fast Gauss transform
R. S. Sampath, H. Sundar and S. Veerapaneni
ACM/IEEE Conference on Supercomputing, 2010. *Best Paper finalist*

We present fast adaptive parallel algorithms to compute the sum of N Gaussians at N points. Direct sequential computation of this sum would take O(N^2) time. The parallel time complexity estimates for our algorithms are O(N/n_p) for uniform point distributions and O(N/n_p log (N/n_p) + n_p log n_p) for nonuniform distributions using n_p CPUs. We incorporate a plane-wave representation of the Gaussian kernel which permits ``diagonal translation''. We use parallel octrees and a new scheme for translating the plane-waves to efficiently handle nonuniform distributions. Computing the transform to six-digit accuracy at 120 billion points took approximately 140 seconds using 4096 cores on the Jaguar supercomputer at the Oak Ridge National Laboratory.

Our implementation is kernel-independent and can handle other ``Gaussian-type'' kernels even when an explicit analytic expression for the kernel is not known. These algorithms form a new class of core computational machinery for solving parabolic PDEs on massively parallel architectures.

The fast generalized Gauss transform
M. Spivak, S. Veerapaneni and L. Greengard
SIAM Journal on Scientific Computing, Vol. 32(5), 2010.

The fast Gauss transform allows for the calculation of the sum of N Gaussians at M points in O(N + M) time. Here, we extend the algorithm to a wider class of kernels, motivated by quadrature issues that arise in using integral equation methods for solving the heat equation on moving domains. In particular, robust high-order product integration methods require convolution with O(q) distinct Gaussian-type kernels in order to obtain qth order accuracy in time. The generalized Gauss transform permits the computation of each of these kernels and, thus, the construction of fast solvers with optimal computational complexity.

We also develop plane-wave representations of these Gaussian-type fields, permitting the ``diagonal translation" version of the Gauss transform to be applied. When the sources and targets lie on a uniform grid, or a hierarchy of uniform grids, we show that the curse of dimensionality (the exponential growth of complexity constants with dimension) can be avoided. Under these conditions, our implementation has a lower operation count than the fast Fourier transform (FFT).

Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method
A. Rahimian, S. Veerapaneni and G. Biros
Journal of Computational Physics, Vol. 229(18), 2010.

...In this paper, we consider two important generalizations: (i) confined flows within arbitrary-shaped stationary/moving geometries and (ii) contrast flows, in which, the fluids enclosed by the vesicles and the ambient fluid can be different from one another. These two problems require solving additional integral
equations and cause modifications to the previous numerical scheme that are nontrivial. Our method does not have severe time-step stability constraints and its computational cost per time step is comparable to that of an explicit scheme. The discretization is pseudo-spectral in space, and multistep BDF in time. We
conduct numerical experiments to investigate the stability, accuracy and the computational cost of the algorithm. Overall, our method achieves several orders of magnitude speed-up compared to standard explicit schemes..

A numerical method for simulating the dynamics of 3D axisymmetric vesicles suspended in viscous flow
S. Veerapaneni, D. Gueyffier, G. Biros and D. Zorin
Journal of Computational Physics, Vol. 228(19), 2009.

We extend ``A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D'', Veerapaneni et al. Journal of Computational Physics, 228(7), 2009 to the case of three dimensional axisymmetric vesicles of spherical or toroidal topology immersed in viscous flows. Although the main components of the algorithm are similar in spirit to the 2D case---spectral approximation in space, semi-implicit time-stepping scheme---the main differences are that the bending and viscous force require new analysis, the linearization for the semi-implicit schemes must be rederived, a fully implicit scheme must be used for the toroidal topology to eliminate a CFL-type restriction, and a novel numerical scheme for the evaluation of the 3D Stokes single-layer potential on an axisymmetric surface is necessary to speed up the calculations. By introducing these novel components, we obtain a time-scheme that experimentally is unconditionally stable, has low cost per time step, and is third-order accurate in time. We present numerical results to analyze the cost and convergence rates of the scheme. To verify the solver, we compare it to a constrained variational approach to compute equilibrium shapes that does not involve interactions with a viscous fluid. To illustrate the applicability of method, we consider a few vesicle-flow interaction problems: the sedimentation of a vesicle, interactions of one and three vesicles with a background Poiseuille flow.

A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D
S. Veerapaneni, D. Gueyffier, D. Zorin and G. Biros
Journal of Computational Physics, Vol. 228(7), 2009.

We present a new method for the evolution of inextensible vesicles immersed in a Stokesian fluid. We use a boundary integral formulation
for the fluid that results in a set of nonlinear integro-differential equations for the vesicle dynamics. The motion of the vesicles is
determined by balancing the nonlocal hydrodynamic forces with the elastic forces due to bending and tension. Numerical simulations of
such vesicle motions are quite challenging. On one hand, explicit time-stepping schemes suffer from a severe stability constraint due to
the stiffness related to high-order spatial derivatives and a milder constraint due to a transport-like stability condition. On the other
hand, an implicit scheme can be expensive because it requires the solution of a set of nonlinear equations at each time step. We
present two semi-implicit schemes that circumvent the severe stability constraints on the time step and whose computational cost per time
step is comparable to that of an explicit scheme. We discretize the equations by using a spectral method in space, and a multistep
third-order accurate scheme in time. We use the fast multipole method (FMM) to efficiently compute vesicle-vesicle interaction forces in a
suspension with a large number of vesicles. We report results from numerical experiments that demonstrate the convergence and algorithmic complexity properties of our scheme.

Analytical and numerical solutions for shapes of quiescent 2D vesicles
S. Veerapaneni, R. Raj, G. Biros and P. K. Purohit
International Journal of Non-linear Mechanics, Vol. 44(3), 2009.

We describe an analytic method for the computation of equilibrium shapes for two-dimensional vesicles characterized by a Helfrich elastic energy. We derive boundary value problems and solve them analytically in terms of elliptic functions and elliptic integrals. We derive solutions by prescribing length and area, or displacements and angle boundary conditions. The solutions are compared to solutions obtained by a boundary integral equation-based numerical scheme. Our method enables the identification of different configurations of deformable vesicles and accurate calculation of their shape, bending moments, tension, and the pressure jump across the vesicle membrane. Furthermore, we perform numerical experiments that indicate that all these configurations are stable minima.

The Chebyshev fast Gauss and nonuniform fast Fourier transforms and their application to the evaluation of distributed heat potentials
S. Veerapaneni and G. Biros
Journal of Computational Physics, Vol. 227 (16), 2008.

We present a method for the fast and accurate computation of distributed heat potentials in two dimensions. The distributed
source is assumed to be given in terms of piecewise space–time Chebyshev polynomials. We discretize uniformly in
time, whereas in space the polynomials are defined on the leaf nodes of a quadtree. The quadtree can vary at each time
step. We combine a product integration rule with fast algorithms (fast heat potentials, nonuniform FFT, fast Gauss transform)
to obtain a high-order accurate method with optimal complexity. If N is the number of time steps, M is the maximum
number of leaf nodes over all the time steps and the input contains a qth-order polynomial representation of f, then,
our method requires O(q^3 MN logM) work to evaluate the heat potential at arbitrary MN space–time target locations. The
overall convergence rate of the method is of order q. We present numerical experiments for q = 4, 8, and 16, and we verify
the theoretical convergence rate of the method. When the solution is sufficiently smooth, the 16th-order variant results in
significant computational savings, even in the case in which we require only a few digits of accuracy.

A fast high-order integral equation solver for the heat equation with moving boundaries in 1D
S. Veerapaneni and G. Biros
SIAM Journal on Scientific Computing, Vol. 29(1), 2007.

We describe a fast high-order accurate method for the solution of the heat equation in domains with moving Dirichlet or Neumann boundaries and distributed forces. We assume that the motion of the boundary is prescribed. Our method extends the work of L. Greengard and J. Strain, "A fast algorithm for the evaluation of heat potentials", Comm. Pure & Applied Math. 1990. Our scheme is based on a time-space Chebyshev pseudo-spectral collocation discretization, which is combined with a recursive product quadrature rule to accurately and efficiently approximate convolutions with the Green's function for the heat equation. We present numerical results that exhibit up to eighth-order convergence rates. Assuming N time steps and M spatial discretization points, the evaluation of the solution of the heat equation at the same number of points in space-time requires O(NM logM) work. Thus, our scheme can be characterized as "fast", that is, it is work-optimal up to a logarithmic factor.