From 16620ff8789b739c0560bcfff1367ecf05fcac88 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 9 Apr 2025 13:16:29 -0500 Subject: [PATCH 1/5] Drop sumpy.kernel.DirectionalTargetDerivative Closes gh-228 --- examples/curve-pot.py | 4 --- sumpy/expansion/multipole.py | 7 ++-- sumpy/kernel.py | 69 ------------------------------------ sumpy/qbx.py | 4 --- sumpy/tools.py | 2 +- 5 files changed, 4 insertions(+), 82 deletions(-) diff --git a/examples/curve-pot.py b/examples/curve-pot.py index c1c113662..c7a4879ab 100644 --- a/examples/curve-pot.py +++ b/examples/curve-pot.py @@ -36,10 +36,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/multipole.py b/sumpy/expansion/multipole.py index ef440cf01..753a98819 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -100,10 +100,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 81556bf92..d9dab1bca 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -79,7 +79,6 @@ .. autoclass:: AxisTargetDerivative .. autoclass:: AxisSourceDerivative .. autoclass:: DirectionalSourceDerivative -.. autoclass:: DirectionalTargetDerivative Transforming kernels -------------------- @@ -1143,74 +1142,6 @@ def __repr__(self): return f"{type(self).__name__}({self.inner_kernel!r}, {self.dir_vec_name})" -class DirectionalTargetDerivative(DirectionalDerivative): - directional_kind = "tgt" - target_array_name = "targets" - - def get_code_transformer(self): - from sumpy.codegen import VectorComponentRewriter - vcr = VectorComponentRewriter([self.dir_vec_name]) - from pymbolic.primitives import Variable - via = _VectorIndexAdder(self.dir_vec_name, (Variable("itgt"),)) - - inner_transform = self.inner_kernel.get_code_transformer() - - def transform(expr): - return via(vcr(inner_transform(expr))) - - return transform - - def postprocess_at_target(self, expr, bvec): - from sumpy.derivative_taker import ( - DifferentiatedExprDerivativeTaker, - diff_derivative_coeff_dict, - ) - from sumpy.symbolic import make_sym_vector as make_sympy_vector - dir_vec = make_sympy_vector(self.dir_vec_name, self.dim) - target_vec = make_sympy_vector(self.target_array_name, self.dim) - - expr = self.inner_kernel.postprocess_at_target(expr, bvec) - - # bvec = tgt - center - if not isinstance(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 += (expr.diff(bvec[axis]) + expr.diff(target_vec[axis])) \ - * dir_vec[axis] - return result - - new_transformation = defaultdict(lambda: 0) - for axis in range(self.dim): - axis_transformation = diff_derivative_coeff_dict( - expr.derivative_coeff_dict, axis, target_vec) - for mi, coeff in axis_transformation.items(): - new_transformation[mi] += coeff * dir_vec[axis] - - return DifferentiatedExprDerivativeTaker(expr.taker, - dict(new_transformation)) - - def get_args(self): - return [ - KernelArgument( - loopy_arg=lp.GlobalArg( - self.dir_vec_name, - None, - shape=(self.dim, "ntargets"), - offset=lp.auto - ), - ), - *self.inner_kernel.get_args() - ] - - def prepare_loopy_kernel(self, loopy_knl): - loopy_knl = self.inner_kernel.prepare_loopy_kernel(loopy_knl) - return lp.tag_array_axes(loopy_knl, self.dir_vec_name, "sep,C") - - mapper_method = "map_directional_target_derivative" - - class DirectionalSourceDerivative(DirectionalDerivative): directional_kind = "src" diff --git a/sumpy/qbx.py b/sumpy/qbx.py index ea3cded78..6491ef7e4 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -512,7 +512,6 @@ def find_jump_term(kernel, arg_provider): AxisTargetDerivative, DerivativeBase, DirectionalSourceDerivative, - DirectionalTargetDerivative, ) tgt_derivatives = [] @@ -522,9 +521,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/tools.py b/sumpy/tools.py index 8251c8bb6..14f76261a 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -298,7 +298,7 @@ def __init__(self, ctx: Any, 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 From 34e1395fe106dfca64cbc650d418dd69f510f6f4 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 9 Apr 2025 13:17:25 -0500 Subject: [PATCH 2/5] Fix derivative taking in line Taylor --- sumpy/expansion/local.py | 9 ++++++--- sumpy/qbx.py | 4 ++-- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 184e0cadd..4697efb89 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 - return sym.Add(*( - coeffs[self.get_storage_index(i)] / math.factorial(i) - for i in self.get_coefficient_identifiers())) + + # 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())) 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/qbx.py b/sumpy/qbx.py index 6491ef7e4..b3016f3a1 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -105,12 +105,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 = [ From d8ecfb7acfb775214bf3531f65486fba6cbe0f5d Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Sat, 28 Jun 2025 14:28:26 -0500 Subject: [PATCH 3/5] add a test for target derivative (starfish) --- test/test_directionaltargetderivative.py | 174 +++++++++++++++++++++++ 1 file changed, 174 insertions(+) create mode 100644 test/test_directionaltargetderivative.py diff --git a/test/test_directionaltargetderivative.py b/test/test_directionaltargetderivative.py new file mode 100644 index 000000000..fb63ff388 --- /dev/null +++ b/test/test_directionaltargetderivative.py @@ -0,0 +1,174 @@ +import pytest +import numpy as np + +import pyopencl as cl +from arraycontext import flatten +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.discretization import Discretization +from meshmode.discretization.poly_element import InterpolatoryQuadratureSimplexGroupFactory +from meshmode.mesh.generation import make_curve_mesh, starfish + +from pytential.qbx import QBXLayerPotentialSource +from pytential import GeometryCollection, bind, sym + +from sumpy.kernel import AxisTargetDerivative, LaplaceKernel +from sumpy.expansion.local import LineTaylorLocalExpansion +from sumpy.qbx import LayerPotentialMatrixGenerator +from pytools.convergence import EOCRecorder + + +def starfish_parametrization(t, n_arms=5, amplitude=0.25): + """Parametrization (x(θ), y(θ)) = r(θ)(cos(θ), sin(θ)), where r(θ) = 1 + amplitude * sin(n_arms * θ).""" + theta = 2 * np.pi * t + + r = 1 + amplitude * np.sin(n_arms * theta) + dr_dt = amplitude * n_arms * 2 * np.pi * np.cos(n_arms * theta) + + x = r * np.cos(theta) + y = r * np.sin(theta) + + dx_dt = dr_dt * np.cos(theta) - r * np.sin(theta) * 2 * np.pi + dy_dt = dr_dt * np.sin(theta) + r * np.cos(theta) * 2 * np.pi + + jacobian_norm = np.sqrt(dx_dt**2 + dy_dt**2) + + tangent_x = dx_dt / jacobian_norm + tangent_y = dy_dt / jacobian_norm + + normal_x = tangent_y + normal_y = -tangent_x + + coords = np.vstack([x, y]) + tangents = np.vstack([tangent_x, tangent_y]) + normals = np.vstack([normal_x, normal_y]) + + return coords, tangents, normals, jacobian_norm + + +@pytest.mark.parametrize("kernel_type", ["laplace"]) +def test_lpot_dx_jump_relation_convergence(kernel_type): + """Test convergence of jump relations for single layer potential derivatives.""" + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + if kernel_type == "laplace": + knl = LaplaceKernel(2) + else: + raise ValueError(f"Unknown kernel type: {kernel_type}") + + qbx_order = 5 + nelements = [100, 150, 200] + target_order = 5 + upsampling_factor = 5 + + ntargets = 20 + np.random.seed(42) + t = np.random.uniform(0, 1, ntargets) + targets_h, _, targets_normals_h, jac = starfish_parametrization( + t, n_arms=5, amplitude=0.25 + ) + targets = actx.from_numpy(targets_h) + + eocrec = EOCRecorder() + + for nelement in nelements: + mesh = make_curve_mesh(starfish, np.linspace(0, 1, nelement + 1), target_order) + pre_density_discr = Discretization( + actx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order) + ) + + qbx = QBXLayerPotentialSource( + pre_density_discr, + upsampling_factor * target_order, + qbx_order, + fmm_order=False + ) + places = GeometryCollection({"qbx": qbx}, auto_where=('qbx')) + + source_discr = places.get_discretization('qbx', sym.QBX_SOURCE_QUAD_STAGE2) + sources_h = actx.to_numpy(flatten(source_discr.nodes(), actx)).reshape(2, -1) + sources = actx.from_numpy(sources_h) + + dofdesc = sym.DOFDescriptor("qbx", sym.QBX_SOURCE_QUAD_STAGE2) + weights_nodes = bind( + places, + sym.weights_and_area_elements(ambient_dim=2, dim=1, dofdesc=dofdesc) + )(actx) + weights_nodes_h = actx.to_numpy(flatten(weights_nodes, actx)) + + expansion_radii_h = jac / (2 * nelement) + centers_in = actx.from_numpy(targets_h - targets_normals_h * expansion_radii_h) + centers_out = actx.from_numpy(targets_h + targets_normals_h * expansion_radii_h) + + mat_gen_x = LayerPotentialMatrixGenerator( + actx.context, + expansion=LineTaylorLocalExpansion(knl, qbx_order), + source_kernels=(knl,), + target_kernels=(AxisTargetDerivative(0, knl),) + ) + + mat_gen_y = LayerPotentialMatrixGenerator( + actx.context, + expansion=LineTaylorLocalExpansion(knl, qbx_order), + source_kernels=(knl,), + target_kernels=(AxisTargetDerivative(1, knl),) + ) + + _, (mat_in_x,) = mat_gen_x( + actx.queue, + targets=targets, + sources=sources, + expansion_radii=expansion_radii_h, + centers=centers_in, + ) + mat_in_x = actx.to_numpy(mat_in_x) + weighted_mat_in_x = mat_in_x * weights_nodes_h[None, :] + + _, (mat_in_y,) = mat_gen_y( + actx.queue, + targets=targets, + sources=sources, + expansion_radii=expansion_radii_h, + centers=centers_in, + ) + mat_in_y = actx.to_numpy(mat_in_y) + weighted_mat_in_y = mat_in_y * weights_nodes_h[None, :] + + _, (mat_out_x,) = mat_gen_x( + actx.queue, + targets=targets, + sources=sources, + expansion_radii=expansion_radii_h, + centers=centers_out, + ) + mat_out_x = actx.to_numpy(mat_out_x) + weighted_mat_out_x = mat_out_x * weights_nodes_h[None, :] + + _, (mat_out_y,) = mat_gen_y( + actx.queue, + targets=targets, + sources=sources, + expansion_radii=expansion_radii_h, + centers=centers_out, + ) + mat_out_y = actx.to_numpy(mat_out_y) + weighted_mat_out_y = mat_out_y * weights_nodes_h[None, :] + + eval_in = (weighted_mat_in_x.sum(axis=1) * targets_normals_h[0] + + weighted_mat_in_y.sum(axis=1) * targets_normals_h[1]) + eval_out = (weighted_mat_out_x.sum(axis=1) * targets_normals_h[0] + + weighted_mat_out_y.sum(axis=1) * targets_normals_h[1]) + + # check jump relation: S'_int - S'_ext = sigma (=1 for constant density) + jump_error = np.abs(eval_in - eval_out - 1) + + h_max = actx.to_numpy(bind(places, sym.h_max(places.ambient_dim))(actx)) + eocrec.add_data_point(h_max, np.max(jump_error)) + + assert eocrec.order_estimate() > qbx_order - 1 + + +if __name__ == "__main__": + pytest.main([__file__]) \ No newline at end of file From e99cc9372e4469d191d3fdb91415910af542414a Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Sat, 28 Jun 2025 15:33:49 -0500 Subject: [PATCH 4/5] matvec -> LayerPotential --- test/test_directionaltargetderivative.py | 79 +++++++++--------------- 1 file changed, 29 insertions(+), 50 deletions(-) diff --git a/test/test_directionaltargetderivative.py b/test/test_directionaltargetderivative.py index fb63ff388..8bb0a1254 100644 --- a/test/test_directionaltargetderivative.py +++ b/test/test_directionaltargetderivative.py @@ -71,8 +71,13 @@ def test_lpot_dx_jump_relation_convergence(kernel_type): ) targets = actx.from_numpy(targets_h) + from sumpy.qbx import LayerPotential + lplot_dx = LayerPotential(actx.context, expansion=LineTaylorLocalExpansion(knl, qbx_order), + target_kernels=(AxisTargetDerivative(0, knl),), source_kernels=(knl,)) + lplot_dy = LayerPotential(actx.context, expansion=LineTaylorLocalExpansion(knl, qbx_order), + target_kernels=(AxisTargetDerivative(1, knl),), source_kernels=(knl,)) eocrec = EOCRecorder() - + for nelement in nelements: mesh = make_curve_mesh(starfish, np.linspace(0, 1, nelement + 1), target_order) pre_density_discr = Discretization( @@ -97,76 +102,50 @@ def test_lpot_dx_jump_relation_convergence(kernel_type): sym.weights_and_area_elements(ambient_dim=2, dim=1, dofdesc=dofdesc) )(actx) weights_nodes_h = actx.to_numpy(flatten(weights_nodes, actx)) - + strengths = (actx.from_numpy(weights_nodes_h),) + expansion_radii_h = jac / (2 * nelement) centers_in = actx.from_numpy(targets_h - targets_normals_h * expansion_radii_h) centers_out = actx.from_numpy(targets_h + targets_normals_h * expansion_radii_h) - - mat_gen_x = LayerPotentialMatrixGenerator( - actx.context, - expansion=LineTaylorLocalExpansion(knl, qbx_order), - source_kernels=(knl,), - target_kernels=(AxisTargetDerivative(0, knl),) - ) - - mat_gen_y = LayerPotentialMatrixGenerator( - actx.context, - expansion=LineTaylorLocalExpansion(knl, qbx_order), - source_kernels=(knl,), - target_kernels=(AxisTargetDerivative(1, knl),) - ) - - _, (mat_in_x,) = mat_gen_x( + + _, (eval_in_dx,) = lplot_dx( actx.queue, - targets=targets, - sources=sources, - expansion_radii=expansion_radii_h, - centers=centers_in, + targets, sources, centers_in, strengths, + expansion_radii=expansion_radii_h ) - mat_in_x = actx.to_numpy(mat_in_x) - weighted_mat_in_x = mat_in_x * weights_nodes_h[None, :] - _, (mat_in_y,) = mat_gen_y( + _, (eval_in_dy,) = lplot_dy( actx.queue, - targets=targets, - sources=sources, - expansion_radii=expansion_radii_h, - centers=centers_in, + targets, sources, centers_in, strengths, + expansion_radii=expansion_radii_h ) - mat_in_y = actx.to_numpy(mat_in_y) - weighted_mat_in_y = mat_in_y * weights_nodes_h[None, :] - _, (mat_out_x,) = mat_gen_x( + _, (eval_out_dx,) = lplot_dx( actx.queue, - targets=targets, - sources=sources, - expansion_radii=expansion_radii_h, - centers=centers_out, + targets, sources, centers_out, strengths, + expansion_radii=expansion_radii_h ) - mat_out_x = actx.to_numpy(mat_out_x) - weighted_mat_out_x = mat_out_x * weights_nodes_h[None, :] - _, (mat_out_y,) = mat_gen_y( + _, (eval_out_dy,) = lplot_dy( actx.queue, - targets=targets, - sources=sources, - expansion_radii=expansion_radii_h, - centers=centers_out, + targets, sources, centers_out, strengths, + expansion_radii=expansion_radii_h ) - mat_out_y = actx.to_numpy(mat_out_y) - weighted_mat_out_y = mat_out_y * weights_nodes_h[None, :] - eval_in = (weighted_mat_in_x.sum(axis=1) * targets_normals_h[0] + - weighted_mat_in_y.sum(axis=1) * targets_normals_h[1]) - eval_out = (weighted_mat_out_x.sum(axis=1) * targets_normals_h[0] + - weighted_mat_out_y.sum(axis=1) * targets_normals_h[1]) + 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 * targets_normals_h[0] + eval_in_dy * targets_normals_h[1] + eval_out = eval_out_dx * targets_normals_h[0] + eval_out_dy * targets_normals_h[1] + # check jump relation: S'_int - S'_ext = sigma (=1 for constant density) jump_error = np.abs(eval_in - eval_out - 1) h_max = actx.to_numpy(bind(places, sym.h_max(places.ambient_dim))(actx)) eocrec.add_data_point(h_max, np.max(jump_error)) - + assert eocrec.order_estimate() > qbx_order - 1 From 6c025ed5e4671b5d018a3bcc4fc097ceeeec5b12 Mon Sep 17 00:00:00 2001 From: Shawn Lin Date: Sat, 28 Jun 2025 16:21:32 -0500 Subject: [PATCH 5/5] comments for the parametrization + Ruff check --- test/test_directionaltargetderivative.py | 127 +++++++++++++---------- 1 file changed, 73 insertions(+), 54 deletions(-) diff --git a/test/test_directionaltargetderivative.py b/test/test_directionaltargetderivative.py index 8bb0a1254..c484da6c1 100644 --- a/test/test_directionaltargetderivative.py +++ b/test/test_directionaltargetderivative.py @@ -1,104 +1,121 @@ -import pytest -import numpy as np +from __future__ import annotations -import pyopencl as cl -from arraycontext import flatten +import numpy as np +import pytest from meshmode.array_context import PyOpenCLArrayContext from meshmode.discretization import Discretization -from meshmode.discretization.poly_element import InterpolatoryQuadratureSimplexGroupFactory +from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory, +) from meshmode.mesh.generation import make_curve_mesh, starfish - -from pytential.qbx import QBXLayerPotentialSource from pytential import GeometryCollection, bind, sym +from pytential.qbx import QBXLayerPotentialSource -from sumpy.kernel import AxisTargetDerivative, LaplaceKernel -from sumpy.expansion.local import LineTaylorLocalExpansion -from sumpy.qbx import LayerPotentialMatrixGenerator +import pyopencl as cl +from arraycontext import flatten from pytools.convergence import EOCRecorder +from sumpy.expansion.local import LineTaylorLocalExpansion +from sumpy.kernel import AxisTargetDerivative, LaplaceKernel + def starfish_parametrization(t, n_arms=5, amplitude=0.25): - """Parametrization (x(θ), y(θ)) = r(θ)(cos(θ), sin(θ)), where r(θ) = 1 + amplitude * sin(n_arms * θ).""" + """ + Parametrization: + (x(θ), y(θ)) = r(θ)(cos(θ), sin(θ)), r(θ) = 1 + amplitude * sin(n_arms * θ). + It is used to compute normal vectors and expansion radius at different + refinement levels for arbitrary boundary targets. + """ + theta = 2 * np.pi * t - + r = 1 + amplitude * np.sin(n_arms * theta) dr_dt = amplitude * n_arms * 2 * np.pi * np.cos(n_arms * theta) - + x = r * np.cos(theta) y = r * np.sin(theta) - + dx_dt = dr_dt * np.cos(theta) - r * np.sin(theta) * 2 * np.pi dy_dt = dr_dt * np.sin(theta) + r * np.cos(theta) * 2 * np.pi - + jacobian_norm = np.sqrt(dx_dt**2 + dy_dt**2) - + tangent_x = dx_dt / jacobian_norm tangent_y = dy_dt / jacobian_norm - - normal_x = tangent_y + + normal_x = tangent_y normal_y = -tangent_x - + coords = np.vstack([x, y]) tangents = np.vstack([tangent_x, tangent_y]) normals = np.vstack([normal_x, normal_y]) - + return coords, tangents, normals, jacobian_norm -@pytest.mark.parametrize("kernel_type", ["laplace"]) +@pytest.mark.parametrize("kernel_type", ["laplace"]) def test_lpot_dx_jump_relation_convergence(kernel_type): """Test convergence of jump relations for single layer potential derivatives.""" - + cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) - + if kernel_type == "laplace": knl = LaplaceKernel(2) else: raise ValueError(f"Unknown kernel type: {kernel_type}") - + qbx_order = 5 nelements = [100, 150, 200] target_order = 5 upsampling_factor = 5 - + ntargets = 20 - np.random.seed(42) - t = np.random.uniform(0, 1, ntargets) + rng = np.random.default_rng(42) + t = rng.uniform(0, 1, ntargets) targets_h, _, targets_normals_h, jac = starfish_parametrization( t, n_arms=5, amplitude=0.25 ) targets = actx.from_numpy(targets_h) from sumpy.qbx import LayerPotential - lplot_dx = LayerPotential(actx.context, expansion=LineTaylorLocalExpansion(knl, qbx_order), - target_kernels=(AxisTargetDerivative(0, knl),), source_kernels=(knl,)) - lplot_dy = LayerPotential(actx.context, expansion=LineTaylorLocalExpansion(knl, qbx_order), - target_kernels=(AxisTargetDerivative(1, knl),), source_kernels=(knl,)) + 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 nelement in nelements: mesh = make_curve_mesh(starfish, np.linspace(0, 1, nelement + 1), target_order) pre_density_discr = Discretization( actx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order) ) - + qbx = QBXLayerPotentialSource( - pre_density_discr, - upsampling_factor * target_order, - qbx_order, + pre_density_discr, + upsampling_factor * target_order, + qbx_order, fmm_order=False ) - places = GeometryCollection({"qbx": qbx}, auto_where=('qbx')) - - source_discr = places.get_discretization('qbx', sym.QBX_SOURCE_QUAD_STAGE2) + places = GeometryCollection({"qbx": qbx}, auto_where=("qbx")) + + source_discr = places.get_discretization("qbx", sym.QBX_SOURCE_QUAD_STAGE2) sources_h = actx.to_numpy(flatten(source_discr.nodes(), actx)).reshape(2, -1) sources = actx.from_numpy(sources_h) - + dofdesc = sym.DOFDescriptor("qbx", sym.QBX_SOURCE_QUAD_STAGE2) weights_nodes = bind( - places, + places, sym.weights_and_area_elements(ambient_dim=2, dim=1, dofdesc=dofdesc) )(actx) weights_nodes_h = actx.to_numpy(flatten(weights_nodes, actx)) @@ -107,47 +124,49 @@ def test_lpot_dx_jump_relation_convergence(kernel_type): expansion_radii_h = jac / (2 * nelement) centers_in = actx.from_numpy(targets_h - targets_normals_h * expansion_radii_h) centers_out = actx.from_numpy(targets_h + targets_normals_h * expansion_radii_h) - + _, (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 * targets_normals_h[0] + eval_in_dy * targets_normals_h[1] - eval_out = eval_out_dx * targets_normals_h[0] + eval_out_dy * targets_normals_h[1] + + eval_in = eval_in_dx * targets_normals_h[0] + \ + eval_in_dy * targets_normals_h[1] + eval_out = eval_out_dx * targets_normals_h[0] + \ + eval_out_dy * targets_normals_h[1] # check jump relation: S'_int - S'_ext = sigma (=1 for constant density) jump_error = np.abs(eval_in - eval_out - 1) - + h_max = actx.to_numpy(bind(places, sym.h_max(places.ambient_dim))(actx)) eocrec.add_data_point(h_max, np.max(jump_error)) - + assert eocrec.order_estimate() > qbx_order - 1 - - + + if __name__ == "__main__": - pytest.main([__file__]) \ No newline at end of file + pytest.main([__file__])