From 2532fcb660c796080ed8e6bfb1aec5f96dbd90b2 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 17 Mar 2025 11:30:25 -0500 Subject: [PATCH 1/4] Fix incorrect scaling of Yukawa kernel Noticed by Shawn Lin. Co-authored-by: Shawn Lin --- sumpy/kernel.py | 4 ++-- sumpy/version.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 12d390ad9..9dfcf14fe 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -676,7 +676,7 @@ def __init__(self, dim: int, yukawa_lambda_name: str = "lam") -> None: expr = var("hankel_1")(0, var("I")*lam*r) scaling_for_K0 = var("pi")/2*var("I") # noqa: N806 - scaling = -1/(2*var("pi")) * scaling_for_K0 + scaling = 1/(2*var("pi")) * scaling_for_K0 elif dim == 3: # NOTE: to get the expression, we do the following and simplify # 1. express K(1/2, lam r) as a modified spherical Bessel function @@ -684,7 +684,7 @@ def __init__(self, dim: int, yukawa_lambda_name: str = "lam") -> None: # 2. or use (AS 10.2.17) directly expr = var("exp")(-lam*r) / r - scaling = -1/(4 * var("pi")**2) + scaling = 1/(4 * var("pi")) else: raise NotImplementedError(f"unsupported dimension: '{dim}'") diff --git a/sumpy/version.py b/sumpy/version.py index 81d1e3c72..8c4f151f9 100644 --- a/sumpy/version.py +++ b/sumpy/version.py @@ -41,4 +41,4 @@ def _parse_version(version: str) -> tuple[tuple[int, ...], str]: VERSION, VERSION_STATUS = _parse_version(VERSION_TEXT) _GIT_REVISION = find_module_git_revision(__file__, n_levels_up=1) -KERNEL_VERSION = (*VERSION, _GIT_REVISION, 0) +KERNEL_VERSION = (*VERSION, _GIT_REVISION, 1) From 99dd45281c20f33f789de01345ade2c03713b4ce Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 24 Mar 2025 16:33:29 -0500 Subject: [PATCH 2/4] Add get_pde_as_diff_op to doc --- sumpy/kernel.py | 1 + 1 file changed, 1 insertion(+) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 9dfcf14fe..60e519d06 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -206,6 +206,7 @@ class Kernel(ABC): .. automethod:: get_global_scaling_const .. automethod:: get_args .. automethod:: get_source_args + .. automethod:: get_pde_as_diff_op """ dim: int From 1ca7cc68a972a4ebe73b5a82bc1de938bbb55e8e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 11 Sep 2025 17:27:45 -0500 Subject: [PATCH 3/4] Add a test for the jump relation to test kernel scaling https://github.com/inducer/sumpy/pull/224 --- sumpy/test/geometries.py | 125 +++++++++++++++++++++++++++++++++++++ sumpy/test/test_kernels.py | 75 +++++++++++++++++++++- 2 files changed, 199 insertions(+), 1 deletion(-) create mode 100644 sumpy/test/geometries.py diff --git a/sumpy/test/geometries.py b/sumpy/test/geometries.py new file mode 100644 index 000000000..7cc829603 --- /dev/null +++ b/sumpy/test/geometries.py @@ -0,0 +1,125 @@ +from __future__ import annotations + + +__copyright__ = "Copyright (C) 2025 University of Illinois Board of Trustees" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from dataclasses import dataclass +from typing import TYPE_CHECKING + +import numpy as np +import numpy.linalg as la + + +if TYPE_CHECKING: + from numpy.typing import NDArray + + +@dataclass(frozen=True) +class Geometry: + nodes: NDArray[np.floating] + normals: NDArray[np.floating] + weights: NDArray[np.floating] + area_elements: NDArray[np.floating] + + +def make_ellipsoid(a: float = 1, b: float = 0.5, npoints: int = 100): + def map_to_curve(t: NDArray[np.floating]): + t = t*(2*np.pi) + + x = a*np.cos(t) + y = b*np.sin(t) + + w = (np.zeros_like(t)+1)/len(t) + + return x, y, w + + from sumpy.test.curve import CurveGrid + + t = np.linspace(0, 1, npoints, endpoint=False) + x, y, weights = map_to_curve(t) + curve = CurveGrid(x, y) + + return Geometry( + nodes=curve.pos, + normals=curve.normal, + weights=weights, + area_elements=curve.speed, + ) + + +def make_torus( + r_major: float = 1, + r_minor: float = 0.5, + n_major: int = 200, + n_minor: int = 200 + ): + u = np.linspace(0.0, 2.0 * np.pi, n_major, endpoint=False) + v = np.linspace(0.0, 2.0 * np.pi, n_minor, endpoint=False) + u, v = np.meshgrid(u, v, copy=False) + nodes = np.stack([ + np.cos(u) * (r_major + r_minor * np.cos(v)), + np.sin(u) * (r_major + r_minor * np.cos(v)), + r_minor * np.sin(v) + ]) + + def diff2d(ary: NDArray[np.floating]): + import scipy.fftpack as fftpack + return np.array([ + fftpack.diff(ary[idx]) + for idx in np.ndindex(ary.shape[:-1]) + ]) + dnodes_du = np.array( + [diff2d(nodes[i].T).T for i in range(3)]) + dnodes_dv = np.array( + [diff2d(nodes[i]) for i in range(3)]) + + normal = -np.cross(dnodes_du, dnodes_dv, axis=0) + area_el = la.norm(normal, axis=0) + normal /= area_el + + weights = np.zeros_like(area_el) + (2*np.pi)**2/n_major/n_minor + + if 0: + import matplotlib.pyplot as plt + plt.figure().add_subplot(projection="3d") + plt.gca().quiver( + nodes[0], nodes[1], nodes[2], + normal[0], normal[1], normal[2], + length=0.1 + ) + plt.show() + 1/0 # noqa: B018 + + geo = Geometry( + nodes=nodes.reshape(3, -1).copy(), + normals=normal.reshape(3, -1).copy(), + area_elements=area_el.reshape(-1).copy(), + weights=weights.reshape(-1).copy(), + ) + + surface_area_ref = 4*np.pi**2*r_major*r_minor + surface_area = geo.area_elements @ geo.weights + surface_area_err = abs(surface_area - surface_area_ref)/surface_area_ref + assert surface_area_err < 1e-14 + + return geo diff --git a/sumpy/test/test_kernels.py b/sumpy/test/test_kernels.py index e34924299..8584dcf6e 100644 --- a/sumpy/test/test_kernels.py +++ b/sumpy/test/test_kernels.py @@ -32,7 +32,11 @@ import pytest import pytools.obj_array as obj_array -from arraycontext import ArrayContextFactory, pytest_generate_tests_for_array_contexts +from arraycontext import ( + ArrayContextFactory, + PyOpenCLArrayContext, + pytest_generate_tests_for_array_contexts, +) from pytools.convergence import PConvergenceVerifier import sumpy.symbolic as sym @@ -41,6 +45,7 @@ from sumpy.expansion.local import ( H2DLocalExpansion, LinearPDEConformingVolumeTaylorLocalExpansion, + LineTaylorLocalExpansion, VolumeTaylorLocalExpansion, ) from sumpy.expansion.m2l import ( @@ -61,7 +66,9 @@ Kernel, LaplaceKernel, StokesletKernel, + YukawaKernel, ) +from sumpy.test.geometries import make_ellipsoid, make_torus logger = logging.getLogger(__name__) @@ -849,6 +856,72 @@ def test_m2m_compressed_error_helmholtz(actx_factory: ArrayContextFactory, dim, # }}} +@pytest.mark.parametrize(("kernel_cls", "kernel_kwargs"), [ + (LaplaceKernel, {}), + (HelmholtzKernel, {"k": 1}), + (YukawaKernel, {"lam": 1}) + ]) +@pytest.mark.parametrize("dim", [2, 3]) +def test_jump( + actx_factory: ArrayContextFactory, + kernel_cls: type[Kernel], + kernel_kwargs: dict[str, float], + dim: int, + order: int = 3, + ): + # x-ref https://github.com/inducer/sumpy/pull/224 + # fails for Yukawa on then-main + + actx = actx_factory() + if not isinstance(actx, PyOpenCLArrayContext): + pytest.skip() + + if dim == 2: + geo = make_ellipsoid(npoints=1600) + qbx_radius = 0.02 + tol = 3e-7 + + elif dim == 3: + geo = make_torus() + qbx_radius = 0.25 + tol = 1e-3 + + else: + raise ValueError(f"unexpected dim: {dim}") + + center_dist = qbx_radius*np.ones(2) + center_sides = np.array([-1, 1], dtype=np.float64) + + targets = geo.nodes[:, 0][:, np.newaxis] + np.zeros(2) + centers = (geo.nodes[:, 0][:, np.newaxis] + + center_sides * center_dist*geo.normals[:, 0][:, np.newaxis]) + + from sumpy.kernel import DirectionalSourceDerivative + kernel = kernel_cls(dim) + dlp_knl = DirectionalSourceDerivative(kernel) + + from sumpy.qbx import LayerPotential + lpot = LayerPotential(actx.context, + expansion=LineTaylorLocalExpansion(kernel, order=order), + source_kernels=(dlp_knl,), + target_kernels=(kernel,), + value_dtypes=np.complex128,) + + _evt, (y,) = lpot(actx.queue, + actx.from_numpy(targets), + actx.from_numpy(geo.nodes), + actx.from_numpy(centers), + [actx.from_numpy(geo.weights * geo.area_elements)], + expansion_radii=actx.from_numpy(center_dist), + src_derivative_dir=actx.from_numpy(geo.normals), + **kernel_kwargs, + ) + + inside, outside = actx.to_numpy(y) + err = abs((inside-outside) - -1) + assert err < tol, err + + # You can test individual routines by typing # $ python test_kernels.py 'test_p2p(_acf, True)' From 154726e073b5d7250a478b1c8ba54265382b2e3e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 11 Sep 2025 17:27:51 -0500 Subject: [PATCH 4/4] Update baseline --- .basedpyright/baseline.json | 298 ++++++++++++++++++++++++++++++++++++ 1 file changed, 298 insertions(+) diff --git a/.basedpyright/baseline.json b/.basedpyright/baseline.json index a2156a4c0..e8caca34f 100644 --- a/.basedpyright/baseline.json +++ b/.basedpyright/baseline.json @@ -42150,6 +42150,264 @@ } } ], + "./sumpy/test/geometries.py": [ + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 20, + "endColumn": 60, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 21, + "endColumn": 41, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 28, + "endColumn": 38, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 28, + "endColumn": 38, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 28, + "endColumn": 36, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 14, + "endColumn": 22, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 4, + "endColumn": 11, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 4, + "endColumn": 10, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 4, + "endColumn": 11, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 28, + "endColumn": 35, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 8, + "endColumn": 18, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 8, + "endColumn": 24, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 16, + "endColumn": 24, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 26, + "endColumn": 34, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 36, + "endColumn": 44, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 16, + "endColumn": 25, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 27, + "endColumn": 36, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 38, + "endColumn": 47, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 8, + "endColumn": 16, + "lineCount": 1 + } + }, + { + "code": "reportUnusedExpression", + "range": { + "startColumn": 8, + "endColumn": 11, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 14, + "endColumn": 27, + "lineCount": 1 + } + }, + { + "code": "reportUnknownMemberType", + "range": { + "startColumn": 14, + "endColumn": 39, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 14, + "endColumn": 41, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 16, + "endColumn": 30, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 16, + "endColumn": 42, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 16, + "endColumn": 44, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 22, + "endColumn": 37, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 22, + "endColumn": 46, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 22, + "endColumn": 48, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 16, + "endColumn": 31, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 16, + "endColumn": 40, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 16, + "endColumn": 42, + "lineCount": 1 + } + } + ], "./sumpy/test/test_codegen.py": [ { "code": "reportUnknownMemberType", @@ -49406,6 +49664,46 @@ "endColumn": 51, "lineCount": 1 } + }, + { + "code": "reportAny", + "range": { + "startColumn": 4, + "endColumn": 10, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 12, + "endColumn": 19, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 36, + "endColumn": 37, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 4, + "endColumn": 7, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 14, + "endColumn": 35, + "lineCount": 1 + } } ], "./sumpy/test/test_matrixgen.py": [