Fix Bug in P2M method (#24)

* Fix bug in p2m method

* Fix mispelled var name in m2l operator

* Fix nodes examined in downward pass

* Fix failure to add multipole expansions for all neighbors in downward pass

fmm/fmm.py

@@ -135,7 +135,8 @@ def downward_pass(self):
         """Downward pass loop."""
 
         # Pre-order traversal of octree
-        for level in range(2, 1 + self.octree.maximum_level):
+        for level in range(1, 1 + self.octree.maximum_level):
+
             for key in self.octree.non_empty_target_nodes_by_level[level]:
                 index = self.octree.target_node_to_index[key]
 
@@ -156,6 +157,7 @@ def downward_pass(self):
             self.local_to_particle(leaf_node_index)
             self.compute_near_field(leaf_node_index)
 
+
     def particle_to_multipole(self, leaf_node_index):
         """Compute particle to multipole interactions in leaf."""
 
@@ -173,42 +175,45 @@ def particle_to_multipole(self, leaf_node_index):
             self.octree.source_leaf_nodes[leaf_node_index]
             ].indices.update(source_indices)
 
-        # Compute key corresponding to this leaf index, and its parent
-        child_key = self.octree.source_leaf_nodes[leaf_node_index]
-        parent_key = hilbert.get_parent(child_key)
+        # Compute key corresponding to this leaf index
+        leaf_key = self.octree.non_empty_source_nodes[leaf_node_index]
 
-        # Compute center of parent box in cartesian coordinates
-        parent_center = hilbert.get_center_from_key(
-            parent_key, self.octree.center, self.octree.radius
+        # Compute center of leaf box in cartesian coordinates
+        leaf_center = hilbert.get_center_from_key(
+            leaf_key, self.octree.center, self.octree.radius
         )
 
         # Compute expansion, and add to source data
         result = p2m(
             kernel_function=self.kernel_function,
             leaf_sources=leaf_sources,
             order=self.order,
-            center=parent_center,
+            center=leaf_center,
             radius=self.octree.radius,
             maximum_level=self.octree.maximum_level
         )
 
-        self.source_data[child_key].expansion = result.density
+        self.source_data[leaf_key].expansion = result.density
 
     def multipole_to_multipole(self, key):
         """Combine children expansions of node into node expansion."""
 
         # Compute center of parent boxes
+
         parent_center = hilbert.get_center_from_key(
             key, self.octree.center, self.octree.radius)
         parent_level = hilbert.get_level(key)
 
+
         for child in hilbert.get_children(key):
+            # Only going through non-empty child nodes
             if self.octree.source_node_to_index[child] != -1:
 
                 # Compute center of child box in cartesian coordinates
                 child_center = hilbert.get_center_from_key(
                     child, self.octree.center, self.octree.radius
                     )
+
                 child_level = hilbert.get_level(child)
 
                 # Updating indices
@@ -266,7 +271,7 @@ def multipole_to_local(self, source_key, target_key):
             source_equivalent_density=source_equivalent_density
         )
 
-        self.target_data[target_key].expansion = result.density
+        self.target_data[target_key].expansion += result.density
 
     def local_to_local(self, key):
         """Translate local expansion of a node to it's children."""
@@ -448,6 +453,17 @@ def __call__(self, x, y):
         raise NotImplementedError
 
 
+class Identity(Kernel):
+    """Identity operation
+    """
+    @staticmethod
+    def kernel_function(x, y):
+        return np.dot(x, y)
+
+    def __call__(self, x, y):
+        return self.kernel_function(x, y)
+
+
 class Laplace(Kernel):
     """Single layer Laplace kernel
     """
@@ -810,7 +826,7 @@ def m2l(kernel_function,
         alpha=2.95
     )
 
-    tgt_check_surface = surface(
+    tgt_downward_check_surface = surface(
         order=order,
         radius=radius,
         level=target_level,
@@ -820,17 +836,20 @@ def m2l(kernel_function,
 
     kernel_tc2te = gram_matrix(
         kernel_function,
-        tgt_check_surface,
+        tgt_downward_check_surface,
         tgt_downward_equivalent_surface
     )
 
     kernel_se2tc = gram_matrix(
-        kernel_function, src_upward_equivalent_surface, tgt_check_surface)
+        kernel_function,
+        src_upward_equivalent_surface,
+        tgt_downward_check_surface
+    )
 
     # Invert gram matrix with SVD
-    kernel_se2te_inv = pseudo_inverse(kernel_tc2te)
+    kernel_tc2te_inv = pseudo_inverse(kernel_tc2te)
 
-    m2l_matrix = np.matmul(kernel_se2te_inv, kernel_se2tc)
+    m2l_matrix = np.matmul(kernel_tc2te_inv, kernel_se2tc)
 
     tgt_equivalent_density = np.matmul(m2l_matrix, source_equivalent_density)
 

fmm/test/test_fmm.py

@@ -101,7 +101,7 @@ def test_p2p(targets, sources, source_densities, expected):
 def test_p2m(n_level_octree, order):
     """Test with single-layer Laplace kernel"""
 
-    octree = n_level_octree(maximum_level=1)
+    octree = n_level_octree(maximum_level=0)
     sources = octree.sources
     source_densities = np.ones(shape=(len(sources)))
 
@@ -110,19 +110,16 @@ def test_p2m(n_level_octree, order):
         order=order,
         center=octree.center,
         radius=octree.radius,
-        maximum_level=1,
+        maximum_level=0,
         leaf_sources=octree.sources
         )
 
-    distant_point = np.array([[10, 0, 0]])
+    distant_point = np.array([[1e4, 0, 0]])
 
     direct = p2p(laplace, distant_point, sources, source_densities)
     equivalent = p2p(laplace, distant_point, result.surface, result.density)
 
-    for i in range(len(equivalent.density)):
-        a = equivalent.density[i]
-        b = direct.density[i]
-        assert np.isclose(a, b, rtol=1e-1)
+    assert np.isclose(equivalent.density, direct.density, rtol=1e-2)
 
     # Test that the surfaces are not equal, as expected
     assert ~np.array_equal(equivalent.surface, direct.surface)
@@ -170,7 +167,7 @@ def test_m2m(n_level_octree, order):
 
     parent_equivalent_surface, parent_equivalent_density = result.surface, result.density
 
-    distant_point = np.array([[1000, 0, 0]])
+    distant_point = np.array([[1e4, 0, 0]])
 
     child_result = p2p(
         laplace, distant_point, child_equivalent_surface, child_equivalent_density
@@ -180,7 +177,7 @@ def test_m2m(n_level_octree, order):
         laplace, distant_point, parent_equivalent_surface, parent_equivalent_density
     )
 
-    assert np.isclose(child_result.density[0], parent_result.density[0], atol=1e-3)
+    assert np.isclose(child_result.density, parent_result.density, rtol=1e-2)
 
     # Test that the surfaces are not equal, as expected
     assert ~np.array_equal(parent_result.surface, child_result.surface)
@@ -357,12 +354,14 @@ def test_upward_pass(level, order, n_level_octree):
         source_densities=np.ones(len(octree.sources))
     )
 
+    print(direct_result)
+
     assert np.isclose(direct_result.density, multipole_result.density, rtol=1.5e-1)
 
 
 def test_downward_pass(n_level_octree, order):
 
-    level = 1
+    level = 2
     octree = n_level_octree(level)
 
     fmm = Fmm(octree, order, laplace)
@@ -375,11 +374,9 @@ def test_downward_pass(n_level_octree, order):
         targets=octree.targets,
         sources=octree.sources,
         source_densities=np.ones(len(octree.sources))
-    )
+    ).density
 
     fmm_p2p = np.array([res.density[0] for res in fmm.result_data])
 
-    print(f"direct p2p: {direct_p2p.density}")
-    print(f"fmm p2p: {fmm_p2p}")
-    print(f"error: {fmm_p2p-direct_p2p.density}")
-    assert False
+    for i in range(len(direct_p2p)):
+        assert np.isclose(direct_p2p[i], fmm_p2p[i], rtol=0.1)