diff --git a/.basedpyright/baseline.json b/.basedpyright/baseline.json index 03a00c5e..57caf352 100644 --- a/.basedpyright/baseline.json +++ b/.basedpyright/baseline.json @@ -17747,15 +17747,15 @@ "code": "reportUnknownArgumentType", "range": { "startColumn": 24, - "endColumn": 60, + "endColumn": 64, "lineCount": 3 } }, { "code": "reportUnknownMemberType", "range": { - "startColumn": 23, - "endColumn": 45, + "startColumn": 27, + "endColumn": 49, "lineCount": 1 } }, @@ -33793,46 +33793,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownMemberType", - "range": { - "startColumn": 21, - "endColumn": 36, - "lineCount": 1 - } - }, - { - "code": "reportUnknownMemberType", - "range": { - "startColumn": 51, - "endColumn": 66, - "lineCount": 1 - } - }, - { - "code": "reportOperatorIssue", - "range": { - "startColumn": 16, - "endColumn": 63, - "lineCount": 1 - } - }, - { - "code": "reportOperatorIssue", - "range": { - "startColumn": 42, - "endColumn": 63, - "lineCount": 1 - } - }, - { - "code": "reportUnknownMemberType", - "range": { - "startColumn": 26, - "endColumn": 38, - "lineCount": 1 - } - }, { "code": "reportOperatorIssue", "range": { @@ -34001,14 +33961,6 @@ "lineCount": 1 } }, - { - "code": "reportAssignmentType", - "range": { - "startColumn": 40, - "endColumn": 66, - "lineCount": 1 - } - }, { "code": "reportUnannotatedClassAttribute", "range": { @@ -38081,14 +38033,6 @@ "lineCount": 1 } }, - { - "code": "reportUnusedParameter", - "range": { - "startColumn": 29, - "endColumn": 33, - "lineCount": 1 - } - }, { "code": "reportUnknownParameterType", "range": { @@ -38161,14 +38105,6 @@ "lineCount": 1 } }, - { - "code": "reportAttributeAccessIssue", - "range": { - "startColumn": 36, - "endColumn": 57, - "lineCount": 1 - } - }, { "code": "reportUnknownMemberType", "range": { @@ -38185,14 +38121,6 @@ "lineCount": 3 } }, - { - "code": "reportUnknownArgumentType", - "range": { - "startColumn": 67, - "endColumn": 68, - "lineCount": 1 - } - }, { "code": "reportUnknownMemberType", "range": { @@ -38225,6 +38153,30 @@ "lineCount": 1 } }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 36, + "endColumn": 43, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 62, + "endColumn": 66, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 68, + "endColumn": 74, + "lineCount": 1 + } + }, { "code": "reportUnknownParameterType", "range": { @@ -39569,30 +39521,6 @@ "lineCount": 1 } }, - { - "code": "reportUnknownMemberType", - "range": { - "startColumn": 12, - "endColumn": 34, - "lineCount": 1 - } - }, - { - "code": "reportUnknownMemberType", - "range": { - "startColumn": 21, - "endColumn": 34, - "lineCount": 1 - } - }, - { - "code": "reportAttributeAccessIssue", - "range": { - "startColumn": 28, - "endColumn": 34, - "lineCount": 1 - } - }, { "code": "reportUnknownArgumentType", "range": { @@ -41719,6 +41647,14 @@ } ], "./sumpy/test/geometries.py": [ + { + "code": "reportAny", + "range": { + "startColumn": 4, + "endColumn": 8, + "lineCount": 1 + } + }, { "code": "reportUnknownArgumentType", "range": { @@ -52176,6 +52112,88 @@ } } ], + "./sumpy/test/test_target_deriv.py": [ + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 35, + "endColumn": 45, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 35, + "endColumn": 45, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 36, + "endColumn": 47, + "lineCount": 1 + } + }, + { + "code": "reportUnknownArgumentType", + "range": { + "startColumn": 36, + "endColumn": 47, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 8, + "endColumn": 15, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 8, + "endColumn": 16, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 8, + "endColumn": 18, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 28, + "endColumn": 50, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 37, + "endColumn": 55, + "lineCount": 1 + } + }, + { + "code": "reportAny", + "range": { + "startColumn": 44, + "endColumn": 54, + "lineCount": 1 + } + } + ], "./sumpy/test/test_tools.py": [ { "code": "reportUnusedImport", diff --git a/examples/curve-pot.py b/examples/curve-pot.py index ae6613e9..10175851 100644 --- a/examples/curve-pot.py +++ b/examples/curve-pot.py @@ -37,10 +37,6 @@ def process_kernel(knl, what_operator): elif what_operator == "D": from sumpy.kernel import DirectionalSourceDerivative source_knl = DirectionalSourceDerivative(knl) - # DirectionalTargetDerivative (temporarily?) removed - # elif what_operator == "S'": - # from sumpy.kernel import DirectionalTargetDerivative - # knl = DirectionalTargetDerivative(knl) else: raise RuntimeError(f"unrecognized operator '{what_operator}'") diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 184e0cad..ebd376f9 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -149,9 +149,12 @@ def coefficients_from_source(self, kernel, avec, bvec, rscale, sac=None): def evaluate(self, tgt_kernel, coeffs, bvec, rscale, sac=None): # no point in heeding rscale here--just ignore it + + # NOTE: We can't meaningfully apply target derivatives here. + # Instead, this is handled in LayerPotentialBase._evaluate. return sym.Add(*( - coeffs[self.get_storage_index(i)] / math.factorial(i) - for i in self.get_coefficient_identifiers())) + coeffs[self.get_storage_index(i)] / math.factorial(i) + for i in self.get_coefficient_identifiers())) def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, dvec, tgt_rscale, sac=None, m2l_translation_classes_dependent_data=None): diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index e230a2af..45f3bced 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -101,10 +101,9 @@ def evaluate(self, kernel, coeffs, bvec, rscale, sac=None): rscale = 1 base_taker = kernel.get_derivative_taker(bvec, rscale, sac) - # Following is a no-op, but AxisTargetDerivative.postprocess_at_target and - # DirectionalTargetDerivative.postprocess_at_target only handles - # DifferentiatedExprDerivativeTaker and sympy expressions, so we need to - # make the taker a DifferentitatedExprDerivativeTaker instance. + # Following is a no-op, but AxisTargetDerivative.postprocess_at_target + # only handles DifferentiatedExprDerivativeTaker and sympy expressions, + # so we need to make the taker a DifferentitatedExprDerivativeTaker instance. base_taker = DifferentiatedExprDerivativeTaker(base_taker, {tuple([0]*self.dim): 1}) taker = kernel.postprocess_at_target(base_taker, bvec) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 60e519d0..9aadf248 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -130,10 +130,6 @@ .. autoclass:: DirectionalSourceDerivative :show-inheritance: :members: mapper_method,directional_kind -.. autoclass:: DirectionalTargetDerivative - :show-inheritance: - :undoc-members: - :members: mapper_method,directional_kind,target_array_name Transforming kernels -------------------- @@ -1351,86 +1347,6 @@ def replace_inner_kernel(self, new_inner_kernel: Kernel) -> Kernel: return type(self)(new_inner_kernel, dir_vec_name=self.dir_vec_name) -class DirectionalTargetDerivative(DirectionalDerivative): - mapper_method: ClassVar[str] = "map_directional_target_derivative" - directional_kind: ClassVar[Literal["src", "tgt"]] = "tgt" - target_array_name: ClassVar[str] = "targets" - - @override - def get_code_transformer(self) -> Callable[[Expression], Expression]: - from sumpy.codegen import VectorComponentRewriter - vcr = VectorComponentRewriter(frozenset([self.dir_vec_name])) - via = _VectorIndexAdder(self.dir_vec_name, (prim.Variable("itgt"),)) - - inner_transform = self.inner_kernel.get_code_transformer() - - def transform(expr: Expression) -> Expression: - return via(vcr(inner_transform(expr))) - - return transform - - @overload - def postprocess_at_target( - self, expr: sym.Expr, bvec: sp.Matrix, - ) -> sym.Expr: ... - - @overload - def postprocess_at_target( - self, expr: ExprDerivativeTaker, bvec: sp.Matrix, - ) -> DifferentiatedExprDerivativeTaker: ... - - @override - def postprocess_at_target( - self, expr: sym.Expr | ExprDerivativeTaker, bvec: sp.Matrix, - ) -> sym.Expr | DifferentiatedExprDerivativeTaker: - dir_vec = sym.make_sym_vector(self.dir_vec_name, self.dim) - target_vec = sym.make_sym_vector(self.target_array_name, self.dim) - - inner_expr = self.inner_kernel.postprocess_at_target(expr, bvec) - - # bvec = tgt - center - if not isinstance(inner_expr, DifferentiatedExprDerivativeTaker): - result = 0 - for axis in range(self.dim): - # Since `bvec` and `tgt` are two different symbolic variables - # need to differentiate by both to get the correct answer - result += ( - (inner_expr.diff(bvec[axis]) + inner_expr.diff(target_vec[axis])) - * dir_vec[axis]) - - assert isinstance(result, sym.Expr) - return result - - new_transformation: DerivativeCoeffDict = defaultdict(lambda: 0) - for axis in range(self.dim): - axis_transformation = diff_derivative_coeff_dict( - inner_expr.derivative_coeff_dict, axis, target_vec) - for mi, coeff in axis_transformation.items(): - new_transformation[mi] += coeff * dir_vec[axis] - - return DifferentiatedExprDerivativeTaker( - inner_expr.taker, dict(new_transformation)) - - @override - def get_args(self) -> Sequence[KernelArgument]: - return [ - KernelArgument( - loopy_arg=lp.GlobalArg( - self.dir_vec_name, - None, - shape=(self.dim, "ntargets"), - offset=lp.auto - ), - ), - *self.inner_kernel.get_args() - ] - - @override - def prepare_loopy_kernel(self, loopy_knl: lp.TranslationUnit) -> lp.TranslationUnit: - loopy_knl = self.inner_kernel.prepare_loopy_kernel(loopy_knl) - return lp.tag_array_axes(loopy_knl, self.dir_vec_name, "sep,C") - - class DirectionalSourceDerivative(DirectionalDerivative): mapper_method: ClassVar[str] = "map_directional_source_derivative" directional_kind: ClassVar[Literal["src", "tgt"]] = "src" @@ -1607,10 +1523,6 @@ def map_axis_source_derivative( self, kernel: AxisSourceDerivative) -> ResultT: return self.rec(kernel.inner_kernel) - def map_directional_target_derivative( - self, kernel: DirectionalTargetDerivative) -> ResultT: - return self.rec(kernel.inner_kernel) - def map_directional_source_derivative( self, kernel: DirectionalSourceDerivative) -> ResultT: return self.rec(kernel.inner_kernel) @@ -1641,11 +1553,6 @@ def map_axis_source_derivative(self, kernel: AxisSourceDerivative) -> Kernel: def map_target_point_multiplier(self, kernel: TargetPointMultiplier) -> Kernel: return type(kernel)(kernel.axis, self.rec(kernel.inner_kernel)) - def map_directional_target_derivative( - self, kernel: DirectionalTargetDerivative) -> Kernel: - return type(kernel)(self.rec(kernel.inner_kernel), - dir_vec_name=kernel.dir_vec_name) - def map_directional_source_derivative( self, kernel: DirectionalSourceDerivative) -> Kernel: return type(kernel)(self.rec(kernel.inner_kernel), @@ -1671,11 +1578,6 @@ def map_axis_target_derivative(self, kernel: AxisTargetDerivative) -> Kernel: class TargetDerivativeRemover(AxisTargetDerivativeRemover): """Removes all target derivatives from the kernel.""" - @override - def map_directional_target_derivative( - self, kernel: DirectionalTargetDerivative) -> Kernel: - return self.rec(kernel.inner_kernel) - class SourceDerivativeRemover(AxisSourceDerivativeRemover): """Removes all source derivatives from the kernel.""" diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 8d6092eb..2864dd69 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -107,12 +107,12 @@ def _expand(self, sac, avec, bvec, rscale, isrc): def _evaluate(self, sac, avec, bvec, rscale, expansion_nr, coefficients): from sumpy.expansion.local import LineTaylorLocalExpansion tgt_knl = self.target_kernels[expansion_nr] - if isinstance(tgt_knl, LineTaylorLocalExpansion): + if isinstance(self.expansion, LineTaylorLocalExpansion): # In LineTaylorLocalExpansion.evaluate, we can't run # postprocess_at_target because the coefficients are assigned # symbols and postprocess with a derivative will make them zero. # Instead run postprocess here before the coefficients are assigned. - coefficients = [tgt_knl.postprocess_at_target(coeff, bvec) for + coefficients = [tgt_knl.postprocess_at_target(coeff, avec) for coeff in coefficients] assigned_coeffs = [ @@ -516,7 +516,6 @@ def find_jump_term(kernel, arg_provider): AxisTargetDerivative, DerivativeBase, DirectionalSourceDerivative, - DirectionalTargetDerivative, ) tgt_derivatives = [] @@ -526,9 +525,6 @@ def find_jump_term(kernel, arg_provider): if isinstance(kernel, AxisTargetDerivative): tgt_derivatives.append(kernel.axis) kernel = kernel.kernel - elif isinstance(kernel, DirectionalTargetDerivative): - tgt_derivatives.append(kernel.dir_vec_name) - kernel = kernel.kernel elif isinstance(kernel, AxisSourceDerivative): src_derivatives.append(kernel.axis) kernel = kernel.kernel diff --git a/sumpy/test/geometries.py b/sumpy/test/geometries.py index 7cc82960..65b243e1 100644 --- a/sumpy/test/geometries.py +++ b/sumpy/test/geometries.py @@ -42,21 +42,15 @@ class Geometry: 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 +def _geometry_from_curve_nodes( + x: NDArray[np.floating], + y: NDArray[np.floating] + ): + npts: int + npts, = x.shape + weights = np.ones(npts, dtype=np.float64)/npts 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( @@ -67,6 +61,25 @@ def map_to_curve(t: NDArray[np.floating]): ) +def make_ellipsoid(a: float = 1, b: float = 0.5, npoints: int = 100): + t = np.linspace(0, 1, npoints, endpoint=False) + x = a*np.cos(2*np.pi*t) + y = b*np.sin(2*np.pi*t) + + return _geometry_from_curve_nodes(x, y) + + +def make_starfish(n_arms: int = 5, amplitude: float = 0.25, npoints: int = 100): + t = np.linspace(0, 1, npoints, endpoint=False) + theta = 2 * np.pi * t + + r = 1 + amplitude * np.sin(n_arms * theta) + + x = r * np.cos(theta) + y = r * np.sin(theta) + return _geometry_from_curve_nodes(x, y) + + def make_torus( r_major: float = 1, r_minor: float = 0.5, diff --git a/sumpy/test/test_target_deriv.py b/sumpy/test/test_target_deriv.py new file mode 100644 index 00000000..cff0b69d --- /dev/null +++ b/sumpy/test/test_target_deriv.py @@ -0,0 +1,153 @@ +from __future__ import annotations + + +__copyright__ = """ +Copyright (C) 2025 Shawn Lin +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. +""" + +import sys +from typing import TYPE_CHECKING + +import numpy as np +import pytest + +from arraycontext import ( + PyOpenCLArrayContext, + pytest_generate_tests_for_array_contexts, +) +from pytools.convergence import EOCRecorder + +from sumpy.array_context import ( # noqa: F401 + PytestPyOpenCLArrayContextFactory, + _acf, # pyright: ignore[reportUnusedImport] +) +from sumpy.expansion.local import LineTaylorLocalExpansion +from sumpy.kernel import AxisTargetDerivative, Kernel, LaplaceKernel +from sumpy.test.geometries import make_starfish + + +if TYPE_CHECKING: + from arraycontext import ArrayContextFactory + +pytest_generate_tests = pytest_generate_tests_for_array_contexts([ + PytestPyOpenCLArrayContextFactory, + ]) + + +@pytest.mark.parametrize("knl", [LaplaceKernel(2)]) +def test_lpot_dx_jump_relation_convergence( + actx_factory: ArrayContextFactory, + knl: Kernel): + """Test convergence of jump relations for single layer potential derivatives.""" + + actx = actx_factory() + if not isinstance(actx, PyOpenCLArrayContext): + pytest.skip() + + qbx_order = 5 + + ntargets = 20 + target_geo = make_starfish(npoints=ntargets) + targets_h = target_geo.nodes + targets = actx.from_numpy(targets_h) + + from sumpy.qbx import LayerPotential + expansion = LineTaylorLocalExpansion(knl, qbx_order) + lplot_dx = LayerPotential( + actx.context, + expansion=expansion, + target_kernels=(AxisTargetDerivative(0, knl),), + source_kernels=(knl,) + ) + lplot_dy = LayerPotential( + actx.context, + expansion=expansion, + target_kernels=(AxisTargetDerivative(1, knl),), + source_kernels=(knl,) + ) + eocrec = EOCRecorder() + + for nsources in [320, 640, 1280, 2560]: + source_geo = make_starfish(npoints=nsources) + sources = actx.from_numpy(source_geo.nodes) + + weights_nodes_h = source_geo.area_elements * source_geo.weights + weights_nodes = actx.from_numpy(weights_nodes_h) + + expansion_radii_h = 4 * target_geo.area_elements / nsources + centers_in = actx.from_numpy( + targets_h - target_geo.normals * expansion_radii_h) + centers_out = actx.from_numpy( + targets_h + target_geo.normals * expansion_radii_h) + + strengths = (weights_nodes,) + _, (eval_in_dx,) = lplot_dx( + actx.queue, + targets, sources, centers_in, strengths, + expansion_radii=expansion_radii_h + ) + + _, (eval_in_dy,) = lplot_dy( + actx.queue, + targets, sources, centers_in, strengths, + expansion_radii=expansion_radii_h + ) + + _, (eval_out_dx,) = lplot_dx( + actx.queue, + targets, sources, centers_out, strengths, + expansion_radii=expansion_radii_h + ) + + _, (eval_out_dy,) = lplot_dy( + actx.queue, + targets, sources, centers_out, strengths, + expansion_radii=expansion_radii_h + ) + + eval_in_dx = actx.to_numpy(eval_in_dx) + eval_in_dy = actx.to_numpy(eval_in_dy) + eval_out_dx = actx.to_numpy(eval_out_dx) + eval_out_dy = actx.to_numpy(eval_out_dy) + + eval_in = eval_in_dx * target_geo.normals[0] + \ + eval_in_dy * target_geo.normals[1] + eval_out = eval_out_dx * target_geo.normals[0] + \ + eval_out_dy * target_geo.normals[1] + + # check jump relation: S'_int - S'_ext = sigma (=1 for constant density) + jump_error = np.abs(eval_in - eval_out - 1) + + h_max = 1/nsources + eocrec.add_data_point(h_max, np.max(jump_error)) + + print(eocrec) + assert eocrec.order_estimate() > qbx_order - 1.5 + + +if __name__ == "__main__": + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + pytest.main([__file__]) diff --git a/sumpy/tools.py b/sumpy/tools.py index 35baf405..2a178eb5 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -320,7 +320,7 @@ def __init__(self, ctx: cl.Context, device: Any | None = None) -> None: """ :arg target_kernels: list of :class:`~sumpy.kernel.Kernel` instances, - with :class:`sumpy.kernel.DirectionalTargetDerivative` as + with :class:`sumpy.kernel.AxisTargetDerivative` as the outermost kernel wrappers, if present. :arg source_kernels: list of :class:`~sumpy.kernel.Kernel` instances with :class:`~sumpy.kernel.DirectionalSourceDerivative` as the