August 2-5, 2010
Boundary integral equations are an efficient approach for the scattering of acoustic and electromagnetic waves from periodic arrays (gratings) of obstacles. The standard way to periodize is to replace the free-space Green's function kernel with its quasi-periodic cousin. This idea has two drawbacks: i) the quasi-periodic Green's function diverges (does not exist) for parameter families known as Wood's anomalies, even though the scattering problem remains well-posed, and ii) the lattice sum representation of the quasi-periodic Green's function converges in a disc, thus becomes unwieldy when obstacles have high aspect ratio.
We bypass both problems by means of a new integral representation that relies on the free-space Green's function alone, and adds auxiliary layer potentials on the boundary of the unit cell strip while expanding the linear system to enforce quasi-periodicity. The result is a (slightly larger) 2nd kind system that achieves spectral accuracy, is immune to Wood's anomalies, avoids lattice sums, handles large aspect ratios, and couples to existing scattering code. By careful summing of neighboring images, obstacles may intersect the unit cell walls.
If time, we will discuss a similar robust new approach to the photonic crystal band structure eigenvalue problem.
Joint work with Leslie Greengard (NYU).
In an approach similar to Ewald’s method for evaluating lattice sums, we split the application of Helmholtz Green’s functions between the spatial and the Fourier domains and, for any finite accuracy, approximate their kernels.
In the spatial domain we use a near optimal linear combination of decaying Gaussians with positive coefficients and, in the Fourier domain, a multiplication by a band-limited kernel obtained by using new quadratures appropriate for the singularity in the Fourier domain. Applying this approach to the free space and the quasi-periodic Green’s functions, as well as those with Dirichlet, Neumann or mixed boundary conditions on simple domains, we obtain fast algorithms in dimensions two and three for computing volumetric integrals involving these Green's functions.
This is a joint work with Chris Kurcz amd Lucas Monzon.
The fast multipole method (FMM) is an efficient algorithm for what is
known as "N-body problems." I will present a new scalable algorithm
and a new implementation of the kernel-independent fast multipole
method, in which both distributed memory parallelism (using MPI) and
shared memory/SIMD parallelism (via GPU acceleration) are employed. I
will conclude my talk by discussing the direct numerical simulation of
blood flow in the Stokes regime using the FMM. I will describe
simulations with 200 million red blood cells, an improvement of four
orders of magnitude over previous results.
George Biros holds Associate Professor appointments with the Schools
of Computational Science and Engineering at Georgia Tech and The
Wallace H. Coulter Department of Biomedical Engineering at Georgia
Tech and Emory University. Prior to joining Georgia Tech, he was an
assistant professor in Mechanical Engineering and Applied Mechanics,
Bioengineering and Computer and Information Science at the University
of Pennsylvania. He received his BS in Mechanical Engineering from
Aristotle University Greece (1995), his MS in Biomedical Engineering
from Carnegie Mellon (1996), and his PhD in Computational Science and
Engineering also from Carnegie Mellon (2000). He was a postdoctoral
associate at the Courant Institute from 2000 to 2003.
We describe a collection of techniques which allow
for the fast and accurate solution of boundary integral
equations on two-dimensional domains whose boundaries
have corner points. Our approach has two key advantages
over existing and recently suggested schemes: (1) it
does not require a prior analytic estimates for solutions,
and (2) many aspects of the scheme generalize readily
to singular three-dimensional domains.
We show how to construct higher-order (10 and above) finite-difference schemes using
Minimum Sobolev Norm techniques that lead to O(N2) condition number discretizations of second-order elliptic PDEs. Numerical experiments comparing with standard FEM solvers show the benefit of the new approach. We then discuss fast numerical methods to convert the discrete equations into discretized integral equations that can have bounded condition number, and hence achieve even higher numerical accuracy.
Sparse solution to an underdetermined linear integral equation is
the central problem for a broad range of applications - scattering,
sensing, imaging, machine learning, signal and image processing, data
analysis and compression, model reduction, optimal control and
design. We will introduce a weak formulation of the problem and
construct its sparse solution by a nonlinear process - the design
of a Gaussian quadrature for the kernel of the integral equation.
We will present a systematic method to solve the resulting quadrature
problem by eigen decomposition.
Solving electromagnetics problem is a challenging task, especially when the structure is multi-scale. These structures are often encountered in circuits in electronic packaging, small antenna designs, as well as small sensor designs. A challenging problem in computational electromagnetics is in solving multi-scale problems in the low frequency regime, especially in the regime between statics and electrodynamics. When the wavelength is much longer than the size of the structure, the physics of the electromagnetic field resembles those of circuits, and hence, this is the circuit physics regime. When the wavelength is sizeable compared to the structure, than wave physics becomes important, and it is important that a simulation method can capture the wave physics interaction. When a structure is multi-scale, and has parts that are small compared to wavelength, but at the same time, is on the order of wavelength, then both circuit physics and wave physics are important. A simulation method has to capture both physics. We will discuss the use of the equivalence principle algorithm (EPA) to capture the multi-scale physics of complex structures. In this method, complex structures are partitioned into parts by the use of equivalence surfaces. The interaction of electromagnetic field with structures within the equivalence surface is done through scattering operators working via the equivalence currents on the equivalence surfaces. The solution within the equivalence surface can be obtained by various numerical methods. Then the interaction between equivalence surfaces is obtained via the use of translation operators. When accelerated with the mixed-form fast multipole method, large multi-scale problems can be solved in this manner. We will also discuss the augmented electric field integral equation (A-EFIE) approach in solving the low-frequency breakdown problem as encountered in circuits in electronic packaging. In this method, the EFIE is augmented with an additional charge unknown, and an additional continuity equation relating the charge to the current. The resultant equation, after proper frequency normalization, is frequency stable down to very low frequency. This belongs to a generalized saddle-point method, and apparently it does not suffer from the low-frequency breakdown, but it does have the low-frequency inaccuracy problem. We will discuss the use of the perturbation method to derive accurate solutions when the low-frequency inaccuracy problem occurs. We will also discuss the hybridization of EPA and A-EFIE to tackle some multi-scale problems.
Joint work with Cris Cecka and Pierre-David Letourneau (Stanford University).
The fast multipole method (FMM) is a technique allowing the fast calculation of sums in O(N) or O(N ln N) steps with some prescribed tolerance. Original FMMs required analytical expansions of the kernel, for example using spherical harmonics or Taylor expansions. In recent years, the range of applicability and the ease of use of FMMs has been considerably extended by the introduction of black box or kernel independent techniques. In these approaches, the user only provides a subroutine to numerically calculate the kernel K. This allows changing the definition of K with practically no change to the computer program. In this talk we will present a novel method that attempts to address some of the shortcomings of existing kernel independent FMMs. This method is applicable to all kernels that are analytic in some appropriate region. An important property of the method is the fact that the multipole-to-local operator is diagonal, that is for p multipole coefficients, the cost of applying the operator is O(p). The computational cost is O(N) for non-oscillatory kernels and O(N ln N) for oscillatory kernels. The approach is based on Cauchy's integral formula. We will present a numerical analysis of the convergence, and methods to find the optimal parameters in the FMM. Numerical results will be presented for some benchmark calculations to demonstrate the accuracy as a function of the number of multipole coefficients, and the computational cost of the different steps.
The butterfly algorithm is a robust alternative to the FFT for
computing certain oscillatory integrals in a fast and accurate manner. In
this approach low-rank interactions are updated in a hierarchical fashion
up and down quadtrees. We review the method, its expected accuracy, and
present an application to synthetic aperture radar imaging. Joint work
with Lexing Ying.
I will discuss recent joint work with Leslie Greengard on a new representation of solutions to the time harmonic Maxwell Equations. Using this representation we reduce the solution of the classical problems of scattering off of a smooth bounded interface to the solution of Fredholm integral equations of second kind on the interface. What distinguishes our representation is that it does not have any spurious "interior" resonances, or suffer from low frequency breakdown. It also reveals some interesting topological features of the time harmonic equations at non-zero frequencies.
The general form of the boundary integral representation formula
for a first order linear system in convervation law form is
The Frechet derivative of this representation is then studied
with the aid
of the parametric derivative formulae for the normal boundary
n ds (resp. n dS or n dV) for boundaries in R2
(resp. R3 or
is shown how the proper handling of the surface variation
recovery of the obvious representation formula the parametric
of the field vector.
Poisson-Boltzmann (PB) equation based implicit solvent model can greatly reduce the computational cost in simulating solvated biomolecular systems by applying the mean field approximation in permittivities and capturing the mobile ions with Boltzmann distribution. However, solving PB equation suffers many numerical difficulties ranging from discontinuous permittivities and electrostatic field across the dielectric interface, the infinite boundary conditions, and charge singularities. In this project, we provide an efficient and accurate numerical algorithm, which adopts a well-conditioned boundary integral equation to handle these difficulties while accelerates the Krylov subspace based iterative methods such as GMRES with treecode. This Cartesian coordinates based treecode is an O(N*log(N)) scheme with properties of easy implementation, efficient memory usage O(N) , and straightforward parallelization. Benchmark testing results on Kirkwood sphere plus simulations of biomoleculars in various sizes are provided to demonstrate the accuracy, efficiency and robustness of the present methods.
In this talk, we present a simple and efficient method for rotating a
spherical harmonic expansion. This is a well-studied problem, arising
in classical scattering theory, quantum mechanics and numerical
analysis, usually addressed through the explicit construction of the
Wigner rotation matrices. Existing fast algorithms, based on
recurrence relations, are subject to a variety of instabilities,
limiting the effectiveness of the approach for expansions of high
We show that rotation can be carried out easily and stably
through "pseudospectral" projection, without ever constructing the
matrix entries themselves. In the simplest version of the method,
projection is carried out on the equator of the rotated sphere. If
only the lowest angular modes are required, the algorithm can be
further accelerated by using a sequence of constant latitude circles.
This is joint work with Leslie Greengard.
We have developed new versions of the FMM in three dimensions for both static and dynamic problems that are reasonably efficient, easy to implement, and fully adaptive. By combining these schemes with suitable quadrature codes, we have constructed "triangle-based" libraries for layer potentials that permit
the rapid deployment of robust integral equation solvers in complex geometry. We will describe some of
the core algorithmic ideas as well as some applications in electrodynamics and acoustics. This is joint work with Zydrunas Gimbutas.
An exact representation is presented for the field inside a sphere
(the observation sphere) due to primary sources enclosed by a second sphere
(the source sphere). The regions bounded by the two spheres have no common
points. The field of the primary sources is expressed in terms of Gaussian
beams whose branch-cut disks are centered in the source sphere. The
expansion coefficients for the standing spherical waves in the observation
sphere are expressed in terms of the output of Gaussian-beam receivers,
whose branch-cut disks are centered in the observation sphere. In this
configuration the patterns of the transmitting and receiving beams
"multiply" to produce a higher directivity than is usually seen with
Gaussian beams. This leads to a fast method for computing matrix-vector
multiplications in scattering calculations, as will be illustrated for a
Dirichlet square plate.