Clarification regarding 3D FMM translation operators

Question

I am implementing an adaptive 3D FMM with the "basic" $O(p^4)$ translation operators. I am looking for clarification on the multipole-to-multipole (M2M) translation operator. I will explain my understanding: In equation 6.26, the arguments $\alpha,\beta$ of the multipole are the coordinates of the child box with respect to the coordinates of the parent box shifted to the origin. Furthermore, the multipole moments $O_{j-n}^{k-m}$ are those of the child, which account for all particles contained in the child and were calculated using (6) where the arguments of the multipole $\alpha_i,\beta_i$ are the coordinates of a particle in the child box with respect to the child box shifted to the origin.

In my implementation, I am testing the M2M operator by evaluating the multipole expansion (eq. 5) at a point sufficiently far from the parent of a child box using (1) the multipole moments of the child box and (2) the multipole moments of the child box translated to the parent box. That is, the only particles under consideration in the parent box are those that are also contained in the child box. My understanding is that the resulting potentials are analytically exact; the only error between these values should be attributed to floating-point error. However, the relative error that my test returns is quite high, see below (I repeat steps (1) and (2) for each of the 4 boxes created in an adaptive refinement).

I wonder if this error is attributed to my understanding of the fundamental theorems as described above? I have tested other aspects of the codebase that I believe could also affect the error.

Answer 1

First, in an FMM implementation, you are always working in some coordinate frame relative to the center of some box, hence the $P-Q$ vector defining the primed coordinate system above. Make sure this is correct in your implementation.

Second, your test is not well-designed for debugging since it looks at an end product rather than the real quantity of interest (multipole coefficients); you should really be doing the following. Calculate the multipole expansion coefficients $M_j^k$ of the parent box directly from the charges contained within. Then, calculate the child multipole expansion coefficients $O_j^k$, apply your M2M operator to obtain $\tilde{M}_j^k$, the parent-level multipole coefficients. These should match $M_j^k$ to many decimal places, if not machine precision. Start with a single-term expansion and verify your implementation gets the constant term correct. Then move to the linear, quadratic, etc. terms to validate your implementation.

Another thing to consider - if your boxes are too small, you will need to apply a scaling for numerical stability. Use boxes with $\mathcal{O}(1)$ diameter for debugging.