35. Fast Multipole Methods

Many simulations in physics require to compute global particle-particle interactions, for example gravity potentials

\[ V(x_j) = - \sum_{i=1}^N \frac{G m_i}{\| x_j - y_i \| } \qquad \forall j, 1 \leq j \leq M, \]

where for \(1 \leq i \leq N\), \(m_i\) are point-masses in \(y_i\). This potential as well it’s gradients must be evaluated in target points \(x_j\). A direct evaluation requires \(O(N M)\) operations. With fast multipole methods one can reduce that to \(O((N+M) \log \varepsilon^{-1})\) operations, where \(\varepsilon\) is the required accuracy.

The FMM - algorithm is considered to be one of the top-10 algorithms of the 20th century: top 10 algorithms

• Leslie Greengard, Vladimir Rokhlin: A fast algorithm for particle simulations, Journal of computational physics 73 (1987)

• Nail A. Gumerov, Ramani Duraiswami: Fast multipole methods for the Helmholtz equation in three dimensions, Elsevier (2004)

We are going to introduce the FMM for the Helmholtz equation, where the summation kernel is Green’s function for the the 3D Helmholtz operator \(\Delta + \kappa^2 I\), with wave number \(\kappa\) $\( k(x,y) = \frac{e^{i \kappa |x-y|} }{4 \pi |x-y|}. \)$

35.1. A quick idea:

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
r = WorkPlane().RectangleC(20,20).Face()
mesh = Mesh(OCCGeometry(r).GenerateMesh(maxh=1))

kappa = 2*pi
def kernel(x,y): 
    dist = Norm(x-y)
    return exp(1j*kappa*dist) / (4*pi*dist)

We start with one charge in the origin, and visualize the radiating field:

pnt1 = CF( (0,0,0) )
pos = CF( (x,y,z) )

Draw (kernel(pnt1, pos), mesh, min=-0.01,max=0.01, animate_complex=True, order=6);

Next we put 3 charges close to the origin, and look at the total potential. The observation suggests to combine the 3 charges to one in the origin:

pnt1 = CF( (0,0,0) )
pnt2 = CF( (0,0.25,0) )
pnt3 = CF( (0,-0.25,0) )
pos = CF( (x,y,z) )
pot = kernel(pnt1,pos) + kernel(pnt2,pos) + kernel(pnt3,pos)
Draw (pot, mesh, min=-0.01,max=0.01, animate_complex=True, order=6);

Two charges of opposite signs form a dipole:

pnt1 = CF( (0,0.25,0) )
pnt2 = CF( (0,-0.25,0) )
pos = CF( (x,y,z) )
pot = kernel(pnt1,pos) - kernel(pnt2,pos) 
Draw (pot, mesh, min=-0.1,max=0.1, animate_complex=True, order=6);

and 4 charges form a quadrupole:

pnt1 = CF( (0.25,0.25,0) )
pnt2 = CF( (0.25,-0.25,0) )
pnt3 = CF( (-0.25,0.25,0) )
pnt4 = CF( (-0.25,-0.25,0) )
pos = CF( (x,y,z) )
pot = kernel(pnt1,pos) - kernel(pnt2,pos) - kernel(pnt3,pos) + kernel(pnt4,pos)
Draw (pot, mesh, min=-0.1,max=0.1, animate_complex=True, order=6);

35.2. A multipole function

NGSolve provides a multipole coefficient-function. This is a sum of point potential, dipole, triplepole, quadrupole etc. up to a specified order. One can add point charges, and it computes the approximation of the potential by a multipole:

from ngsolve.fem import SingularMultiPoleCF
from ngsolve.bla import Vec3D

mp = SingularMultiPoleCF(10, kappa, Vec3D(0,0,0))
mp.AddCharge (Vec3D(0,0.25,0), 1)
mp.AddCharge (Vec3D(0,-0.25,0), -1)

Draw (mp, mesh, min=-0.1,max=0.1, animate_complex=True, order=6);

This approximation cannot work everywhere: The original potential is singular at the given charge, but the multipole approximation is a smooth function. Now let’s look closer to the sources. We observer a divergent multipole sum inside a ball, which radius is given by the maximal distance of the charge to the origin of the multipole. Outside the radius we observe a convergent series, where the convergence is faster if we are further from the origin.

mp = SingularMultiPoleCF(20, kappa/5, Vec3D(0,0,0))
mp.AddCharge (Vec3D(0,1,0), 1)
mp.AddCharge (Vec3D(0,-1,0), -1)

Draw (mp, mesh, min=-0.1,max=0.1, animate_complex=True, order=3);

35.3. Regular expansion

Away from the origin of the multipole, the function is smooth. We introduce regular approximations on balls contained in the domain of convergence of the singular multipole. We call them regular multipoles. The Transform method of a multipole object computes the expansion coefficients of the target multipole.

• We have a charge at (2,0,0)

• We approximate its potential by a singular multipole centered at (0,0,0). Domain of convergence is outside a ball of radius 2

• The we approximate that singular multipole by a regular multipole centered at (5,0,0). It converges inside a ball of radius 3 (centered at (5,0,0,)) to the potential of the given charge. Outside that ball we observe divergence of the regular multipole.

from ngsolve.fem import RegularMultiPoleCF

mp = SingularMultiPoleCF(20, kappa, Vec3D(0,0,0))
mp.AddCharge (Vec3D(2,0,0), 3)

regmp = RegularMultiPoleCF(40, kappa, Vec3D(5,0,0) )
mp.Transform(regmp)

Draw (mp, mesh, min=-0.1,max=0.1, animate_complex=True, order=5);
Draw (regmp, mesh, min=-0.1,max=0.1, animate_complex=True, order=5);

We see that we cannot approximate the potentials everywhere by sindular or regular multipoles. But, we can fill the space by balls (or complements of balls), and compute multipole coefficients for each ball, and use that multipole to represent the function.