From 7ed95392ceff5c89bb28145f7753e6c25ef4808f Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 21 Oct 2019 13:12:26 -0500 Subject: [PATCH 1/9] Revert "Fix biharmonic kernel expressions" This reverts commit 7cf3197b8378fa16739f3c740bd80e9ead04765d. --- sumpy/kernel.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 3936c518..3edef3c6 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -447,14 +447,13 @@ class BiharmonicKernel(ExpressionKernel): init_arg_names = ("dim",) def __init__(self, dim=None): - # See https://arxiv.org/abs/1202.1811 r = pymbolic_real_norm_2(make_sym_vector("d", dim)) if dim == 2: - expr = r**2 * (var("log")(r) - 1) - scaling = -1/(8*var("pi")) + expr = r**2 * var("log")(r) + scaling = 1/(8*var("pi")) elif dim == 3: expr = r - scaling = 1/(8*var("pi")) + scaling = 1 # FIXME: Unknown else: raise RuntimeError("unsupported dimensionality") -- GitLab From 624289956778a4cdd3de27c7f6468a23e038c959 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Sun, 17 Nov 2019 20:34:04 -0600 Subject: [PATCH 2/9] Support multiple source derivatives in multipole expansion --- sumpy/expansion/multipole.py | 45 +++++++++++++++++++++++++----------- 1 file changed, 31 insertions(+), 14 deletions(-) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 47e33a8c..7bcd0ed5 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -47,6 +47,17 @@ class MultipoleExpansionBase(ExpansionBase): pass +def cartesian_product(*args): + if len(args) == 1: + for arg in args[0]: + yield (arg,) + return + first = args[:-1] + for prod in cartesian_product(*first): + for i in args[-1]: + yield prod + (i,) + + # {{{ volume taylor class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): @@ -64,27 +75,33 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): rscale = 1 if isinstance(kernel, DirectionalSourceDerivative): - if kernel.get_base_kernel() is not kernel.inner_kernel: - raise NotImplementedError("more than one source derivative " - "not supported at present") - from sumpy.symbolic import make_sym_vector - dir_vec = make_sym_vector(kernel.dir_vec_name, kernel.dim) + + dir_vecs = [] + tmp_kernel = kernel + while isinstance(tmp_kernel, DirectionalSourceDerivative): + dir_vecs.append(make_sym_vector(tmp_kernel.dir_vec_name, kernel.dim)) + tmp_kernel = tmp_kernel.inner_kernel + + if kernel.get_base_kernel() is not tmp_kernel: + raise NotImplementedError("Unknown kernel wrapper.") + + nderivs = len(dir_vecs) coeff_identifiers = self.get_full_coefficient_identifiers() result = [0] * len(coeff_identifiers) - for idim in range(kernel.dim): - for i, mi in enumerate(coeff_identifiers): - if mi[idim] == 0: + for i, mi in enumerate(coeff_identifiers): + for value in cartesian_product(*[range(kernel.dim)]*nderivs): + prod = 1 + derivative_mi = list(mi) + for nderivative, deriv_dim in enumerate(value): + prod *= -derivative_mi[deriv_dim]*dir_vecs[nderivative][deriv_dim] + derivative_mi[deriv_dim] -= 1 + if any(v < 0 for v in derivative_mi): continue + result[i] += mi_power(avec, derivative_mi) * prod - derivative_mi = tuple(mi_i - 1 if iaxis == idim else mi_i - for iaxis, mi_i in enumerate(mi)) - - result[i] += ( - - mi_power(avec, derivative_mi) * mi[idim] - * dir_vec[idim]) for i, mi in enumerate(coeff_identifiers): result[i] /= (mi_factorial(mi) * rscale ** sum(mi)) else: -- GitLab From 5f0916cfdfb25379b5f070a4d2e2d2ad221f222f Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Sun, 17 Nov 2019 22:09:11 -0600 Subject: [PATCH 3/9] Use cartesian_product from pytools --- requirements.txt | 1 + sumpy/expansion/multipole.py | 12 +----------- 2 files changed, 2 insertions(+), 11 deletions(-) diff --git a/requirements.txt b/requirements.txt index 8e48bc1c..8bd95658 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,3 +6,4 @@ git+https://github.com/inducer/pyopencl git+https://gitlab.tiker.net/inducer/boxtree git+https://github.com/inducer/loopy git+https://github.com/inducer/pyfmmlib +git+https://github.com/inducer/pytools diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 7bcd0ed5..4f6ffe3e 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -29,6 +29,7 @@ from sumpy.symbolic import vector_xreplace from sumpy.expansion import ( ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, HelmholtzConformingVolumeTaylorExpansion) +from pytools import cartesian_product import logging logger = logging.getLogger(__name__) @@ -47,17 +48,6 @@ class MultipoleExpansionBase(ExpansionBase): pass -def cartesian_product(*args): - if len(args) == 1: - for arg in args[0]: - yield (arg,) - return - first = args[:-1] - for prod in cartesian_product(*first): - for i in args[-1]: - yield prod + (i,) - - # {{{ volume taylor class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): -- GitLab From b4dee08e4409f79b3bb91b674d22c01c91930c26 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Sun, 17 Nov 2019 23:47:33 -0600 Subject: [PATCH 4/9] Fix tests --- requirements.txt | 2 +- sumpy/expansion/multipole.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/requirements.txt b/requirements.txt index 8bd95658..efe91f94 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,9 +1,9 @@ numpy sympy==1.1.1 +git+https://github.com/inducer/pytools git+https://github.com/inducer/pymbolic git+https://github.com/inducer/islpy git+https://github.com/inducer/pyopencl git+https://gitlab.tiker.net/inducer/boxtree git+https://github.com/inducer/loopy git+https://github.com/inducer/pyfmmlib -git+https://github.com/inducer/pytools diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 4f6ffe3e..7b08778c 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -86,7 +86,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): prod = 1 derivative_mi = list(mi) for nderivative, deriv_dim in enumerate(value): - prod *= -derivative_mi[deriv_dim]*dir_vecs[nderivative][deriv_dim] + prod *= -derivative_mi[deriv_dim] + prod *= dir_vecs[nderivative][deriv_dim] derivative_mi[deriv_dim] -= 1 if any(v < 0 for v in derivative_mi): continue -- GitLab From 8dd56294cd31f696b58a1bf1d3745cdd4230a425 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 18 Nov 2019 01:30:45 -0600 Subject: [PATCH 5/9] Update pytools on conda env --- .test-conda-env-py3.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index e0434a18..aa8ae80a 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -18,6 +18,7 @@ dependencies: - pip - pip: + - git+https://github.com/inducer/pytools - git+https://gitlab.tiker.net/inducer/boxtree - git+https://github.com/inducer/pymbolic - git+https://github.com/inducer/loopy -- GitLab From d9c2017b0a76184a552d245bb4b83b2a1db10598 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 25 Nov 2019 18:31:42 -0600 Subject: [PATCH 6/9] Add explanations and warnings --- sumpy/expansion/multipole.py | 9 +++++++-- sumpy/kernel.py | 6 ++++++ 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 7b08778c..2d6cbf40 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -82,10 +82,15 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): result = [0] * len(coeff_identifiers) for i, mi in enumerate(coeff_identifiers): - for value in cartesian_product(*[range(kernel.dim)]*nderivs): + # One source derivative is the dot product of the gradient and + # directional vector. + # For multiple derivatives, gradient of gradients is taken. + # For eg: in 3D, 2 source derivatives gives 9 terms and + # cartesian_product below enumerates these 9 terms. + for deriv_terms in cartesian_product(*[range(kernel.dim)]*nderivs): prod = 1 derivative_mi = list(mi) - for nderivative, deriv_dim in enumerate(value): + for nderivative, deriv_dim in enumerate(deriv_terms): prod *= -derivative_mi[deriv_dim] prod *= dir_vecs[nderivative][deriv_dim] derivative_mi[deriv_dim] -= 1 diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 3edef3c6..cfed294a 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -449,11 +449,17 @@ class BiharmonicKernel(ExpressionKernel): def __init__(self, dim=None): r = pymbolic_real_norm_2(make_sym_vector("d", dim)) if dim == 2: + # Ref: Farkas, Peter. Mathematical foundations for fast algorithms + # for the biharmonic equation. Diss. University of Chicago, + # Dept. of Mathematics, 1989. expr = r**2 * var("log")(r) scaling = 1/(8*var("pi")) elif dim == 3: expr = r scaling = 1 # FIXME: Unknown + from warnings import warn + warn("scaling factor for Biharmonic 3D is unknown.", + stacklevel=2) else: raise RuntimeError("unsupported dimensionality") -- GitLab From f1beacdfcb38f30e186cdeec12394caee3f72f2c Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 25 Nov 2019 18:31:54 -0600 Subject: [PATCH 7/9] Add test --- test/test_kernels.py | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/test/test_kernels.py b/test/test_kernels.py index 3fcb195b..ce5160e5 100644 --- a/test/test_kernels.py +++ b/test/test_kernels.py @@ -37,13 +37,15 @@ from sumpy.expansion.multipole import ( VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion, VolumeTaylorMultipoleExpansionBase, LaplaceConformingVolumeTaylorMultipoleExpansion, - HelmholtzConformingVolumeTaylorMultipoleExpansion) + HelmholtzConformingVolumeTaylorMultipoleExpansion, + BiharmonicConformingVolumeTaylorMultipoleExpansion) from sumpy.expansion.local import ( VolumeTaylorLocalExpansion, H2DLocalExpansion, LaplaceConformingVolumeTaylorLocalExpansion, - HelmholtzConformingVolumeTaylorLocalExpansion) + HelmholtzConformingVolumeTaylorLocalExpansion, + BiharmonicConformingVolumeTaylorLocalExpansion) from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative, - DirectionalSourceDerivative) + DirectionalSourceDerivative, BiharmonicKernel) from pytools.convergence import PConvergenceVerifier import logging @@ -121,6 +123,15 @@ def test_p2p(ctx_getter, exclude_self): (HelmholtzKernel(2), H2DLocalExpansion), (HelmholtzKernel(2), H2DMultipoleExpansion), + (DirectionalSourceDerivative(BiharmonicKernel(2), "dir_vec"), + VolumeTaylorMultipoleExpansion), + (DirectionalSourceDerivative(BiharmonicKernel(2), "dir_vec"), + VolumeTaylorLocalExpansion), + (DirectionalSourceDerivative(BiharmonicKernel(2), "dir_vec"), + BiharmonicConformingVolumeTaylorMultipoleExpansion), + (DirectionalSourceDerivative(BiharmonicKernel(2), "dir_vec"), + BiharmonicConformingVolumeTaylorLocalExpansion), + (HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorMultipoleExpansion), (HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorLocalExpansion), (HelmholtzKernel(2, allow_evanescent=True), @@ -193,7 +204,7 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative): box_source_counts_nonchild = np.array([nsources], dtype=np.int32) extra_source_kwargs = extra_kwargs.copy() - if with_source_derivative: + if isinstance(knl, DirectionalSourceDerivative): alpha = np.linspace(0, 2*np.pi, nsources, np.float64) dir_vec = np.vstack([np.cos(alpha), np.sin(alpha)]) extra_source_kwargs["dir_vec"] = dir_vec @@ -322,6 +333,10 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative): slack += 1 grad_slack += 2 + if isinstance(base_knl, DirectionalSourceDerivative): + slack += 1 + grad_slack += 2 + if isinstance(base_knl, HelmholtzKernel): if base_knl.allow_evanescent: slack += 0.5 -- GitLab From f1a1566882e54f3d1301972729782a5e9a133c7f Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 25 Nov 2019 20:34:33 -0600 Subject: [PATCH 8/9] Fix citation --- sumpy/kernel.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index cfed294a..b0e245e7 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -450,8 +450,8 @@ class BiharmonicKernel(ExpressionKernel): r = pymbolic_real_norm_2(make_sym_vector("d", dim)) if dim == 2: # Ref: Farkas, Peter. Mathematical foundations for fast algorithms - # for the biharmonic equation. Diss. University of Chicago, - # Dept. of Mathematics, 1989. + # for the biharmonic equation. Technical Report 765, + # Department of Computer Science, Yale University, 1990. expr = r**2 * var("log")(r) scaling = 1/(8*var("pi")) elif dim == 3: -- GitLab From 1cec3c4c3bc1fe3b13de0a6ec7f4155eed609c42 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 25 Nov 2019 21:38:25 -0600 Subject: [PATCH 9/9] Fix testing --- test/test_kernels.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/test/test_kernels.py b/test/test_kernels.py index ce5160e5..834d8233 100644 --- a/test/test_kernels.py +++ b/test/test_kernels.py @@ -37,13 +37,11 @@ from sumpy.expansion.multipole import ( VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion, VolumeTaylorMultipoleExpansionBase, LaplaceConformingVolumeTaylorMultipoleExpansion, - HelmholtzConformingVolumeTaylorMultipoleExpansion, - BiharmonicConformingVolumeTaylorMultipoleExpansion) + HelmholtzConformingVolumeTaylorMultipoleExpansion) from sumpy.expansion.local import ( VolumeTaylorLocalExpansion, H2DLocalExpansion, LaplaceConformingVolumeTaylorLocalExpansion, - HelmholtzConformingVolumeTaylorLocalExpansion, - BiharmonicConformingVolumeTaylorLocalExpansion) + HelmholtzConformingVolumeTaylorLocalExpansion) from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative, DirectionalSourceDerivative, BiharmonicKernel) from pytools.convergence import PConvergenceVerifier @@ -127,10 +125,6 @@ def test_p2p(ctx_getter, exclude_self): VolumeTaylorMultipoleExpansion), (DirectionalSourceDerivative(BiharmonicKernel(2), "dir_vec"), VolumeTaylorLocalExpansion), - (DirectionalSourceDerivative(BiharmonicKernel(2), "dir_vec"), - BiharmonicConformingVolumeTaylorMultipoleExpansion), - (DirectionalSourceDerivative(BiharmonicKernel(2), "dir_vec"), - BiharmonicConformingVolumeTaylorLocalExpansion), (HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorMultipoleExpansion), (HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorLocalExpansion), -- GitLab