Fast boundary element method for the Laplace and Helmholtz equations in three dimensions

The boundary element method (BEM) is a well known technique for solution of boundary value problems appearing in various fields of physics, such as fluid and solid mechanics, electromagnetism, acoustics, etc. To use this method one should discretize the boundary of the domain with M elements (or N~M vertices) and solve a discrete analog of boundary integral equation to determine potential and its normal derivative on the boundary. Knowledge of these quantities allows one to construct the field at any spatial point.

Despite many advantages of the BEM, which enables handling of objects of complex shape, it has some computational drawbacks. One of them is substantial memory (O(N²)) required to store the matrices, which include all pairwise interactions between the elements. Another is the speed, which depends on the number of operations to solve the system. Conventionally, to solve a system of N linear equations one should spend O(N³) operations. Substantial improvement here can be achieved using iterative methods. Using standard matrix-vector multiplication which is of computational cost O(N²) one can solve with some accuracy the same system of linear equations using only O(N_iterN²) operations, where N_iter is the number of iterations, which is presumably much smaller than N. Even this improvement appears to be not sufficient to solve large problems where N is of order tens of thousands or more. The breakthrough method which allows us to solve million size problems on a desktop PC is the Fast Multipole Method (FMM). This method requires much less memory (just O(N)), as it does not need to store all pairwise interactions, and has a computational cost for single matrix-vector product O(N), asit employs factorization and hierarchical data structure.

In the example shown above we put randomly oriented ellipsoids into a cubic mesh and solved a boundary value problem for the Laplace equation. The number of ellipsoids varied from 1 to 1728. The latter configuration had 1,679,616 surface elements and 843,264 vertices (this is not the limit and the numbers are selected just for nice fit to a graph on the right, which has scale up to 1,000,000 vertices). Solution of this large problem took only 28 min 19 s on my desktop PC (3.2 GHz Intel Xeon, 3.5 GB RAM) and required 31 GMRES iterations (accuracy of the solution was controlled) with 48 s per matrix-vector multiplication. As one can see from the performance graph conventional BEM does not fit the memory of my computer for problems of size larger than N=10,000. Moreover these methods are slow. Indeed, assume that there is no memory limitation. In this case we can extrapolate the CPU time required by the direct solution method and the iterative method with matrix storage to find that the first method should take more than 21 years and the second one more than 2 days (which one can afford if 160 Terabytes = 160,000 GB of RAM is available). The performance of the FMM is slightly larger than O(N) (in fact O(N^1.1), because in this example the number of iterations increases with N (not substantially).

We also developed the BEM/FMM technique for the Helmholtz equation, which is the wave equation in the frequency domain. Particularly we were interested with solution of the scattering problems from biological objects, which usually have very complex shape. The problems we solved had up to 2 millions of elements and 1 million of vertices (which fitted to my desktop CPU).

One of the problems is shown below. Here the bunny mesh was exposed to the plane wave and the scattering problem for a sound hard object was solved. A movie showing acoustic pressure was generated based on this solution. This shows the total pressure on the bunny with 32 frames per period of acoustic wave.

Mesh: 132,072 surface elements, 65539 vertices	250 Hz (kD=0.96)
	2.5 kHz (kD=9.6)
	25 kHz (kD=96)

Another problem we have a long story to deal with is multiple scattering problem (see also below). One of computational cases (plane wave scattering from 512 ellipsoids, 249,856 vertices/497,664 elements, kD=144.45) is shown below (this animates acoustic pressure on the ellipsoid surfaces).

©2006, Nail Gumerov

Computation of acoustic scattering from N spheres using multipole reexpansion method

Numerous practical problems of acoustic and electromagnetic wave propagation require computation of the field scattered by multiple objects. Examples include scattering of acoustic waves by objects (e.g., the scattering of sound by humans and the environment), light scattering by clouds and environment, electromagnetic waves in composite materials and the human body, pressure waves in disperse systems (aerosols, emulsions, bubbly liquids), etc.

In many cases the scatterers are spheres, or can be modeled as such. Such modeling is convenient for parametrization of large problems, since each sphere can be characterized by a few quantities such as the coordinates of its center, its radius, and its impedance. This impedance will in general be a complex quantity, and characterizes the absorbing/reflecting properties of the body/surface. For example, we are exploring the modeling of the human head and body using two spheres representing respectively the head and the torso. In fluid mechanical problems, bubbles, droplets, or dust particles can be assumed spherical.

We are interested in computing the solution of multiple scattering from N spheres, with specified impedance boundary conditions at their surfaces.

We developed a method (variant of the T-matrix method) and software that implements it (Journal of the Acoustical Society of America, 112(6), 2002, 2688-2701; and also see technical report). The method is fast enough, since it uses a recursive computation of multipole reexpansion coefficients (Chew, 1992; Gumerov & Duraiswami, 2001).

We compared the results of the computations with numerical and analytical solutions, and demonstrate the computational efficiency of our method in comparison with commercial BEM software. The results show that the developed method compares favorably in both accuracy and speed up of computations (in some cases by several orders of magnitude).

Currently we are working on even faster methods to solve the same problem, which includes implementation of the iterartive solvers (Reflection method and GMRES), and fast matrix-vector multiplication engines (via fast multipole translation operators and fast multipole method). Some preliminary results show that we move in a right direction.

Incident Field Total Field Scattered Field

These pictures demonstrate an example computation with our code MultisphereHelmholtz of scattering of a monochromatic wave from 4 spheres, which centers are arranged in tetrahedron. These spheres are exposed to an incident plane wave that comes from the direction shown by the arrow (the angle of the arrow with the normal to the imaging plane shown in blue is 45 degrees and it belongs to a symmetry plane). In computations we selected parameter ka=15.2 (where k is the wavenumber, k=2p/l). The picture on the left hand side shows distribution of the amplitude of the acoustic field over the sphere surfaces. The imaging plane surves only for visualization of the field. By tracking the lines of constant phase, we can animate the wave propagation in time. This is illustrated by the three shadow pictures on the right hand side. In the absence of spheres the plane wave is imaged as moving shadow strips. The presence of spheres results in diffraction patterns (the middle shadow picture). The difference between the total field and the incident field amplitudes yields the patterns of the scattered field.

ka = 1.6 ka = 2.8 ka = 4.8

This set of pictures illustrate an example computation with MultisphereHelmholtz of scattering of a monochromatic wave from 100 spheres, which centers are randomly uniformly distributed inside a cube and radii also are random numbers (between 0.5 and 1)(see the picture on the left hand side). The incident field is a plane wave which direction is indicated by a small arrow on top of the left picture. This direction is perpendicular to the imaging plane shown by blue dots. Computations were performed for different wavelengths and produced shadow diffraction patterns that are shown on the pictures. It is seen that the scale of the chaotic structures is determined by the wavelength.

If your computer is OK to run 1 MB .gif movie you can click here to view color scattering patterns for this problem. The movie models polychromatic wave scattering (superposition of the three monochromatic waves shown above). Each wavelength is imaged by its own color. We put the above movies ka=1.6 to the red, ka=2.8 to the green and ka=4.8 to the blue channel of the RGB-color representation.

The capabilities of the direct computations, however, are limited, and the problems of size 100 or so, perhaps are the maximum that one can solve on a personal computer. The Fast Multipole Method provides dramatical speed up of computations and enables solution of problems with thousands of spheres on the same PC with the same accuracy. The paper describing this method is published in Journal of the Acoustical Society of America, March, 2005; also see technical report). With this method we computed problems for scattering off more than 100,000 spheres (this is not the limit). To view 3MB. gif movie showing the speckle patterns for 10000 spheres (monochromatic wave ka=0.75, randomly uniform distribution of equal spheres inside a cube, volume fraction 0.2), please click here.

©2003,2005, Nail Gumerov.