Code:realpole

This page contains code related to the JCP article, "Real-space quadrature: A convenient, efficient representation for multipole expansions"

To use the code, you will also need a set of spherical quadrature rules. These were originally from [1], but are included in the package below for ease of use. The code can all be downloaded from [2]:

git clone https://github.com/frobnitzem/realpole

It contains several python modules:

integral.py

Integral.py defines the Lebedev class, instantiated by providing a polynomial order, e.g. <source lang="python"> L = Lebedev(11) </source> it holds quadrature points and weights for polynomials up to (and including) the degree specified. Note that all odd polynomials integrate to zero, so Lebedev(10) = Lebedev(11). The class can integrate a function of a vector by evaluating the function and summing. L.integrate(lambda x: x[0]*x[0]) will numerically integrate <math>x^2</math> over the sphere.

Note that this class needs the quadrature rules in the rules/ subdirectory. Here is its rule loading code: <source lang="python"> cwd = os.path.dirname(os.path.abspath(__file__))

def read_rule(K):

   r = read_matrix("%s/rules/lebedev_%03d.txt"%(cwd,K))

</source>

Integral also contains the e2poly function, which generates Legendre polynomials.

Tests

The test() function in Integral.py runs through the set of all monomials and tests the numerical against the analytical integral. Be sure to use this if you change the quadrature rules.

Poisson.py

The module Poisson.py provides the Poisson class, responsible for working with moment expansions. Notice that its first task is to load up the quadrature rules for polynomials up to order 2*(K-1)+1: <source lang="python"> self.L = Lebedev(2*(K-1)+1) </source> This is because it represents simultaneously all spherical harmonics of order 0, 1, ..., K-1. Integrating products therefore requires rules capable of integrating twice the order of the highest polynomial. The +1 is because all the rules are odd.

At the heart of this class is the moment-generating method, <source lang="python">

   # Calculate the polynomials in series:
   # sum_k q_k L_n(x_k, x_i) [d=1]
   # sum_k q_k L_n(x_i, x_k) [d=-1]
   def mom(self, q, x, d=1): # Hi Mom!
       if d == -1:
           for P in e2poly(self.L.x, x, self.K, -1):
               yield dot(P, q)
       else:
           for P in e2poly(x, self.L.x, self.K, -1):
               yield dot(q, P)

</source>

This is used to find the effective surface charges (weights) in the high-level interface, <source lang="python">

   # If d = 1 (forward),
   # return {w} such that sum q_a*L_n(x_a,y) = sum w_i*L_n(x_i,y)
   # for all K moments L_0, ..., L_{K-1}
   #
   # If d = -1 (reverse),
   # return {w} such that sum q_a*L_n(x,y_a) = sum w_i*L_n(x,y_i)
   # for all K moments L_0, ..., L_{K-1}
   #
   # -- scales as O(K * len(x) * N)
   def solve_moments(self, q,x,R,d=1):

</source>

The code also contains ishift and oshift functions implementing the inner and outer shift operations. See the function description (and usage in mom_t.py) for their documentation.

Tests

The function test(R) will generate a random distribution of point charges and run the following tests:

• rep_err Check that the moment fitting problem was solved by comparing quadrature-fitted vs. exact moments. This test should give a numerical zero.
• pot_err Check the error of representing point charges by their quadrature representation. This test should give decreasing error as the expansion order increases.
• outer_shift Check the error of the outer -> outer shift operation.
• inner_shift Check the error of the inner -> inner shift operation.
• io_shift Check the error of the outer -> inner shift operation.

These last three tests should give error that decreases as the expansion order increases.

The files mom.py and mom_t.py were used to generate the error plots (Figures 2 and 3) in the paper. These are just expanded versions of the tests above.

fmm.py

This is a very basic implementation of the fast multipole method (not used in the paper). It needs some reorganization to make it simpler to experiment with options that tune the error and scaling. At this point, the testing is preliminary and has an error dominated by the size of the interaction list. I would welcome code contributions or suggestions for integrating with PyFMM here!