From 96df67d561da14b3c17e3dfbcbfd88fb7122ff14 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 19 Jul 2019 12:39:39 -0500 Subject: [PATCH 1/5] Improve FMMLIB numerical stability for Laplace 3D --- boxtree/pyfmmlib_integration.py | 3 + test/test_fmm.py | 166 +++++++++++++++++++++++++------- 2 files changed, 135 insertions(+), 34 deletions(-) diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py index 2ca27e3..9a1f225 100644 --- a/boxtree/pyfmmlib_integration.py +++ b/boxtree/pyfmmlib_integration.py @@ -199,6 +199,9 @@ class FMMLibExpansionWrangler(object): # }}} def level_to_rscale(self, level): + if self.dim == 3 and self.eqn_letter == "l": + return 1 + result = self.tree.root_extent * 2 ** -level * self.rscale_factor if abs(result) > 1: result = 1 diff --git a/test/test_fmm.py b/test/test_fmm.py index e673ce9..4b28f0c 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -47,6 +47,52 @@ logger = logging.getLogger(__name__) warnings.simplefilter("ignore", FMMLibRotationDataNotSuppliedWarning) +# {{{ ref fmmlib pot computation + +def get_fmmlib_ref_pot(wrangler, weights, sources_host, targets_host, + helmholtz_k, dipole_vec=None): + dims = sources_host.shape[0] + eqn_letter = "h" if helmholtz_k else "l" + use_dipoles = dipole_vec is not None + + import pyfmmlib + fmmlib_routine = getattr( + pyfmmlib, + "%spot%s%ddall%s_vec" % ( + eqn_letter, + "fld" if dims == 3 else "grad", + dims, + "_dp" if use_dipoles else "")) + + kwargs = {} + if dims == 3: + kwargs["iffld"] = False + else: + kwargs["ifgrad"] = False + kwargs["ifhess"] = False + + if use_dipoles: + if helmholtz_k == 0 and dims == 2: + kwargs["dipstr"] = ( + -weights # pylint:disable=invalid-unary-operand-type + * (dipole_vec[0] + 1j * dipole_vec[1])) + else: + kwargs["dipstr"] = weights + kwargs["dipvec"] = dipole_vec + else: + kwargs["charge"] = weights + if helmholtz_k: + kwargs["zk"] = helmholtz_k + + return wrangler.finalize_potentials( + fmmlib_routine( + sources=sources_host, targets=targets_host, + **kwargs)[0] + ) + +# }}} + + # {{{ fmm interaction completeness test class ConstantOneExpansionWranglerWithFilteredTargetsInTreeOrder( @@ -419,40 +465,8 @@ def test_pyfmmlib_fmm(ctx_factory, dims, use_dipoles, helmholtz_k): logger.info("computing direct (reference) result") - import pyfmmlib - fmmlib_routine = getattr( - pyfmmlib, - "%spot%s%ddall%s_vec" % ( - wrangler.eqn_letter, - "fld" if dims == 3 else "grad", - dims, - "_dp" if use_dipoles else "")) - - kwargs = {} - if dims == 3: - kwargs["iffld"] = False - else: - kwargs["ifgrad"] = False - kwargs["ifhess"] = False - - if use_dipoles: - if helmholtz_k == 0 and dims == 2: - kwargs["dipstr"] = ( - -weights # pylint:disable=invalid-unary-operand-type - * (dipole_vec[0] + 1j * dipole_vec[1])) - else: - kwargs["dipstr"] = weights - kwargs["dipvec"] = dipole_vec - else: - kwargs["charge"] = weights - if helmholtz_k: - kwargs["zk"] = helmholtz_k - - ref_pot = wrangler.finalize_potentials( - fmmlib_routine( - sources=sources_host.T, targets=targets_host.T, - **kwargs)[0] - ) + ref_pot = get_fmmlib_ref_pot(wrangler, weights, sources_host.T, + targets_host.T, helmholtz_k, dipole_vec) rel_err = la.norm(pot - ref_pot, np.inf) / la.norm(ref_pot, np.inf) logger.info("relative l2 error vs fmmlib direct: %g" % rel_err) @@ -508,6 +522,90 @@ def test_pyfmmlib_fmm(ctx_factory, dims, use_dipoles, helmholtz_k): # }}} +# {{{ test fmmlib numerical stability + +@pytest.mark.parametrize("dims", [2, 3]) +@pytest.mark.parametrize("helmholtz_k", [0, 2]) +@pytest.mark.parametrize("order", [35]) +def test_pyfmmlib_numerical_stability(ctx_factory, dims, helmholtz_k, order): + logging.basicConfig(level=logging.INFO) + + from pytest import importorskip + importorskip("pyfmmlib") + + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + nsources = 20 + dtype = np.float64 + + # The input particles are arranged with geometrically decreasing spacing + # along a line, to build a deep tree that stress-tests the translations. + zero = np.zeros(nsources, dtype=dtype) + sources = np.vstack([ + np.array([2**-i for i in range(nsources)], dtype=dtype), + zero, + zero])[:dims] + + targets = sources * (1 + 3 ** -nsources) + + from boxtree import TreeBuilder + tb = TreeBuilder(ctx) + + tree, _ = tb(queue, sources, targets=targets, + max_particles_in_box=2, debug=True) + + assert tree.nlevels >= 15 + + from boxtree.traversal import FMMTraversalBuilder + tbuild = FMMTraversalBuilder(ctx) + trav, _ = tbuild(queue, tree, debug=True) + + trav = trav.get(queue=queue) + + from pyopencl.clrandom import PhiloxGenerator + rng = PhiloxGenerator(queue.context, seed=20) + + weights = rng.uniform(queue, nsources, dtype=np.float64).get() + + from boxtree.pyfmmlib_integration import ( + FMMLibExpansionWrangler, FMMLibRotationData) + + def fmm_level_to_nterms(tree, lev): + return order + + wrangler = FMMLibExpansionWrangler( + trav.tree, helmholtz_k, + fmm_level_to_nterms=fmm_level_to_nterms, + rotation_data=FMMLibRotationData(queue, trav)) + + from boxtree.fmm import drive_fmm + + pot = drive_fmm(trav, wrangler, weights) + assert not np.isnan(pot).any() + + # {{{ ref fmmlib computation + + logger.info("computing direct (reference) result") + + ref_pot = get_fmmlib_ref_pot(wrangler, weights, sources, targets, + helmholtz_k) + + rel_err = la.norm(pot - ref_pot, np.inf) / la.norm(ref_pot, np.inf) + logger.info("relative l2 error vs fmmlib direct: %g" % rel_err) + + if dims == 2: + error_bound = (1/2) ** (1 + order) + else: + error_bound = (3/4) ** (1 + order) + + assert rel_err < error_bound, rel_err + + # }}} + +# }}} + + # {{{ test particle count thresholding in traversal generation @pytest.mark.parametrize("enable_extents", [True, False]) -- GitLab From 5374d9e3c26900843d99c8ca07dbf8d6144e56d6 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 19 Jul 2019 14:30:00 -0500 Subject: [PATCH 2/5] Use a more complicated tree to ensure List 2 is being tested --- test/test_fmm.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/test/test_fmm.py b/test/test_fmm.py index 4b28f0c..af29649 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -536,18 +536,23 @@ def test_pyfmmlib_numerical_stability(ctx_factory, dims, helmholtz_k, order): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - nsources = 20 + nsources = 30 dtype = np.float64 - # The input particles are arranged with geometrically decreasing spacing - # along a line, to build a deep tree that stress-tests the translations. + # The input particles are arranged with geometrically increasing/decreasing + # spacing along a line, to build a deep tree that stress-tests the + # translations. + particle_line = np.array([2**-i for i in range(nsources//2)], dtype=dtype) + particle_line = np.hstack([particle_line, 3 - particle_line]) zero = np.zeros(nsources, dtype=dtype) + sources = np.vstack([ - np.array([2**-i for i in range(nsources)], dtype=dtype), + particle_line, zero, zero])[:dims] - targets = sources * (1 + 3 ** -nsources) + targets = sources.copy() + targets[0] += 5 from boxtree import TreeBuilder tb = TreeBuilder(ctx) @@ -562,11 +567,7 @@ def test_pyfmmlib_numerical_stability(ctx_factory, dims, helmholtz_k, order): trav, _ = tbuild(queue, tree, debug=True) trav = trav.get(queue=queue) - - from pyopencl.clrandom import PhiloxGenerator - rng = PhiloxGenerator(queue.context, seed=20) - - weights = rng.uniform(queue, nsources, dtype=np.float64).get() + weights = np.ones_like(sources[0]) from boxtree.pyfmmlib_integration import ( FMMLibExpansionWrangler, FMMLibRotationData) -- GitLab From 51fb3850c9aae1a6b4a375e881c276de3b8b93a7 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 19 Jul 2019 14:38:00 -0500 Subject: [PATCH 3/5] Move targets closer to sources --- test/test_fmm.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_fmm.py b/test/test_fmm.py index af29649..2bd3f21 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -551,8 +551,7 @@ def test_pyfmmlib_numerical_stability(ctx_factory, dims, helmholtz_k, order): zero, zero])[:dims] - targets = sources.copy() - targets[0] += 5 + targets = sources * (1 + 1e-3) from boxtree import TreeBuilder tb = TreeBuilder(ctx) -- GitLab From 4f30c19748fe6d97f7973b157d326a6c8cddd945 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Tue, 23 Jul 2019 01:25:37 +0200 Subject: [PATCH 4/5] Set rscale to 1/boxsize for L3D --- boxtree/pyfmmlib_integration.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py index 9a1f225..23ae43b 100644 --- a/boxtree/pyfmmlib_integration.py +++ b/boxtree/pyfmmlib_integration.py @@ -199,12 +199,11 @@ class FMMLibExpansionWrangler(object): # }}} def level_to_rscale(self, level): - if self.dim == 3 and self.eqn_letter == "l": - return 1 - result = self.tree.root_extent * 2 ** -level * self.rscale_factor if abs(result) > 1: result = 1 + if self.dim == 3 and self.eqn_letter == "l": + result = 1 / result return result @memoize_method -- GitLab From a4cf92c98b764d8d7ba759f26f2614d4993e8046 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Tue, 23 Jul 2019 01:38:20 +0200 Subject: [PATCH 5/5] Add a comment explaining the rscale flip --- boxtree/pyfmmlib_integration.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py index 23ae43b..937605c 100644 --- a/boxtree/pyfmmlib_integration.py +++ b/boxtree/pyfmmlib_integration.py @@ -203,6 +203,9 @@ class FMMLibExpansionWrangler(object): if abs(result) > 1: result = 1 if self.dim == 3 and self.eqn_letter == "l": + # Laplace 3D uses the opposite convention compared to + # all other cases. + # https://gitlab.tiker.net/inducer/boxtree/merge_requests/81 result = 1 / result return result -- GitLab