diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index e0434a1803073bf6dfa78879021d597e3fdee55c..aa8ae80adf6740131933528c34ecb98d8d7fb0c3 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 diff --git a/requirements.txt b/requirements.txt index 8e48bc1c8b0da2b4d81682839c688c0ac69637d5..efe91f9487cb17f063b456973aabd909f153aae7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,6 @@ 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 diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 47e33a8c76ae8ad5f85e72efcab7415fbd648253..2d6cbf409b73b83023e75731738033cc219cddca 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__) @@ -64,27 +65,39 @@ 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): + # 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(deriv_terms): + 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 + 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: diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 3936c518246aa8f14e54dba2f0b33d5361da04c6..b0e245e7f3987b67d59ab7a3ab7a7c56ac985f03 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -447,14 +447,19 @@ 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")) + # Ref: Farkas, Peter. Mathematical foundations for fast algorithms + # 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: expr = r - scaling = 1/(8*var("pi")) + scaling = 1 # FIXME: Unknown + from warnings import warn + warn("scaling factor for Biharmonic 3D is unknown.", + stacklevel=2) else: raise RuntimeError("unsupported dimensionality") diff --git a/test/test_kernels.py b/test/test_kernels.py index 3fcb195bba0be9cb0b67a4822a45dec616d370aa..834d823321bec0c05b3253e4c8718ff9145b343f 100644 --- a/test/test_kernels.py +++ b/test/test_kernels.py @@ -43,7 +43,7 @@ from sumpy.expansion.local import ( LaplaceConformingVolumeTaylorLocalExpansion, HelmholtzConformingVolumeTaylorLocalExpansion) from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative, - DirectionalSourceDerivative) + DirectionalSourceDerivative, BiharmonicKernel) from pytools.convergence import PConvergenceVerifier import logging @@ -121,6 +121,11 @@ 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), + (HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorMultipoleExpansion), (HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorLocalExpansion), (HelmholtzKernel(2, allow_evanescent=True), @@ -193,7 +198,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 +327,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