diff --git a/fmm/fmm.py b/fmm/fmm.py index a6081aa..af2a556 100644 --- a/fmm/fmm.py +++ b/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,13 +175,12 @@ 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 @@ -187,28 +188,32 @@ def particle_to_multipole(self, leaf_node_index): 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) diff --git a/fmm/test/test_fmm.py b/fmm/test/test_fmm.py index 4282da9..6f04ca2 100644 --- a/fmm/test/test_fmm.py +++ b/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 \ No newline at end of file + for i in range(len(direct_p2p)): + assert np.isclose(direct_p2p[i], fmm_p2p[i], rtol=0.1) \ No newline at end of file