From 93fad2c5f5ec619b24b3dfe0c7519840348eb1af Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 18 Jan 2018 11:22:42 -0600 Subject: [PATCH 1/2] Move unregularized tests to their own file. --- test/test_layer_pot.py | 91 --------------------- test/test_unregularized.py | 162 +++++++++++++++++++++++++++++++++++++ 2 files changed, 162 insertions(+), 91 deletions(-) create mode 100644 test/test_unregularized.py diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index cb8669e4..f6212dd7 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -225,97 +225,6 @@ def test_off_surface_eval_vs_direct(ctx_getter, do_plot=False): # }}} -# {{{ unregularized tests - - -def test_unregularized_with_ones_kernel(ctx_getter): - cl_ctx = ctx_getter() - queue = cl.CommandQueue(cl_ctx) - - nelements = 10 - order = 8 - - mesh = make_curve_mesh(partial(ellipse, 1), - np.linspace(0, 1, nelements+1), - order) - - from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory - - discr = Discretization(cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(order)) - - from pytential.unregularized import UnregularizedLayerPotentialSource - lpot_src = UnregularizedLayerPotentialSource(discr) - - from sumpy.kernel import one_kernel_2d - - expr = sym.IntG(one_kernel_2d, sym.var("sigma"), qbx_forced_limit=None) - - from pytential.target import PointsTarget - op_self = bind(lpot_src, expr) - op_nonself = bind((lpot_src, PointsTarget(np.zeros((2, 1), dtype=float))), expr) - - with cl.CommandQueue(cl_ctx) as queue: - sigma = cl.array.zeros(queue, discr.nnodes, dtype=float) - sigma.fill(1) - sigma.finish() - - result_self = op_self(queue, sigma=sigma) - result_nonself = op_nonself(queue, sigma=sigma) - - assert np.allclose(result_self.get(), 2 * np.pi) - assert np.allclose(result_nonself.get(), 2 * np.pi) - - -def test_unregularized_off_surface_fmm_vs_direct(ctx_getter): - cl_ctx = ctx_getter() - queue = cl.CommandQueue(cl_ctx) - - nelements = 300 - target_order = 8 - fmm_order = 4 - - mesh = make_curve_mesh(WobblyCircle.random(8, seed=30), - np.linspace(0, 1, nelements+1), - target_order) - - from pytential.unregularized import UnregularizedLayerPotentialSource - from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory - - density_discr = Discretization( - cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) - direct = UnregularizedLayerPotentialSource( - density_discr, - fmm_order=False, - ) - fmm = direct.copy( - fmm_level_to_order=lambda kernel, kernel_args, tree, level: fmm_order) - - sigma = density_discr.zeros(queue) + 1 - - fplot = FieldPlotter(np.zeros(2), extent=5, npoints=100) - from pytential.target import PointsTarget - ptarget = PointsTarget(fplot.points) - from sumpy.kernel import LaplaceKernel - - op = sym.D(LaplaceKernel(2), sym.var("sigma"), qbx_forced_limit=None) - - direct_fld_in_vol = bind((direct, ptarget), op)(queue, sigma=sigma) - fmm_fld_in_vol = bind((fmm, ptarget), op)(queue, sigma=sigma) - - err = cl.clmath.fabs(fmm_fld_in_vol - direct_fld_in_vol) - - linf_err = cl.array.max(err).get() - print("l_inf error:", linf_err) - assert linf_err < 5e-3 - -# }}} - - # {{{ test performance data gathering def test_perf_data_gathering(ctx_getter, n_arms=5): diff --git a/test/test_unregularized.py b/test/test_unregularized.py new file mode 100644 index 00000000..c104acd4 --- /dev/null +++ b/test/test_unregularized.py @@ -0,0 +1,162 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2017, 2018 Matt Wala" + +__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 numpy as np +import numpy.linalg as la # noqa +import pyopencl as cl +import pyopencl.clmath # noqa +import pytest # noqa +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +from functools import partial +from meshmode.mesh.generation import ( # noqa + ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle, + make_curve_mesh, NArmedStarfish) +from sumpy.visualization import FieldPlotter +from pytential import bind, sym, norm # noqa + +import logging +logger = logging.getLogger(__name__) + +circle = partial(ellipse, 1) + + +# {{{ ones kernel test + +def test_ones_kernel(ctx_getter): + cl_ctx = ctx_getter() + queue = cl.CommandQueue(cl_ctx) + + nelements = 10 + order = 8 + + mesh = make_curve_mesh(partial(ellipse, 1), + np.linspace(0, 1, nelements+1), + order) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + discr = Discretization(cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(order)) + + from pytential.unregularized import UnregularizedLayerPotentialSource + lpot_src = UnregularizedLayerPotentialSource(discr) + + from sumpy.kernel import one_kernel_2d + + expr = sym.IntG(one_kernel_2d, sym.var("sigma"), qbx_forced_limit=None) + + from pytential.target import PointsTarget + op_self = bind(lpot_src, expr) + op_nonself = bind((lpot_src, PointsTarget(np.zeros((2, 1), dtype=float))), expr) + + with cl.CommandQueue(cl_ctx) as queue: + sigma = cl.array.zeros(queue, discr.nnodes, dtype=float) + sigma.fill(1) + sigma.finish() + + result_self = op_self(queue, sigma=sigma) + result_nonself = op_nonself(queue, sigma=sigma) + + assert np.allclose(result_self.get(), 2 * np.pi) + assert np.allclose(result_nonself.get(), 2 * np.pi) + +# }}} + + +# {{{ on surface direct + +# }}} + + +# {{{ on surface fmm + +# }}} + + +# {{{ off surface fmm versus direct + +def test_off_surface_fmm_vs_direct(ctx_getter): + cl_ctx = ctx_getter() + queue = cl.CommandQueue(cl_ctx) + + nelements = 300 + target_order = 8 + fmm_order = 4 + + mesh = make_curve_mesh(WobblyCircle.random(8, seed=30), + np.linspace(0, 1, nelements+1), + target_order) + + from pytential.unregularized import UnregularizedLayerPotentialSource + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + density_discr = Discretization( + cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) + direct = UnregularizedLayerPotentialSource( + density_discr, + fmm_order=False, + ) + fmm = direct.copy( + fmm_level_to_order=lambda kernel, kernel_args, tree, level: fmm_order) + + sigma = density_discr.zeros(queue) + 1 + + fplot = FieldPlotter(np.zeros(2), extent=5, npoints=100) + from pytential.target import PointsTarget + ptarget = PointsTarget(fplot.points) + from sumpy.kernel import LaplaceKernel + + op = sym.D(LaplaceKernel(2), sym.var("sigma"), qbx_forced_limit=None) + + direct_fld_in_vol = bind((direct, ptarget), op)(queue, sigma=sigma) + fmm_fld_in_vol = bind((fmm, ptarget), op)(queue, sigma=sigma) + + err = cl.clmath.fabs(fmm_fld_in_vol - direct_fld_in_vol) + + linf_err = cl.array.max(err).get() + print("l_inf error:", linf_err) + assert linf_err < 5e-3 + +# }}} + + +# You can test individual routines by typing +# $ python test_layer_pot.py 'test_routine()' + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + import py.test + py.test.cmdline.main([__file__]) + +# vim: fdm=marker -- GitLab From aef116a53ea7e093ff1ba543898414da1749ae4f Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 29 Jan 2018 21:01:26 -0600 Subject: [PATCH 2/2] [WIP] [ci skip] More work towards unregularized self eval with the trapezoid rule. --- pytential/qbx/__init__.py | 10 +++++ pytential/source.py | 42 ++++++++++++------- pytential/symbolic/compiler.py | 29 +++++++++++-- pytential/symbolic/mappers.py | 5 ++- pytential/symbolic/primitives.py | 17 ++++++-- pytential/unregularized.py | 58 ++++++++++++++++++++------ test/test_unregularized.py | 71 ++++++++++++++++++++++++++++---- 7 files changed, 190 insertions(+), 42 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index e05a8063..5698733b 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -56,6 +56,16 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): .. attribute :: qbx_order .. attribute :: fmm_order + .. rubric:: Discretizations + + Supports refinement and upsampling. + + .. attribute:: density_discr + .. attribute:: stage2_density_discr + .. attribute:: quad_stage2_density_discr + .. attribute:: resampler + .. method:: with_refinement + See :ref:`qbxguts` for some information on the inner workings of this. """ diff --git a/pytential/source.py b/pytential/source.py index 6420494e..a740e427 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -94,7 +94,7 @@ class PointPotentialSource(PotentialSource): return result @memoize_method - def get_p2p(self, kernels): + def get_p2p(self, kernels, exclude_self=False): # needs to be separate method for caching from pytools import any @@ -104,8 +104,9 @@ class PointPotentialSource(PotentialSource): value_dtype = self.real_dtype from sumpy.p2p import P2P - p2p = P2P(self.cl_context, - kernels, exclude_self=False, value_dtypes=value_dtype) + p2p = P2P( + self.cl_context, kernels, exclude_self=exclude_self, + value_dtypes=value_dtype) return p2p @@ -151,16 +152,7 @@ class PointPotentialSource(PotentialSource): # {{{ layer potential source class LayerPotentialSourceBase(PotentialSource): - """A discretization of a layer potential using panel-based geometry, with - support for refinement and upsampling. - - .. rubric:: Discretizations - - .. attribute:: density_discr - .. attribute:: stage2_density_discr - .. attribute:: quad_stage2_density_discr - .. attribute:: resampler - .. method:: with_refinement + """A discretization of a layer potential using panel-based geometry. .. rubric:: Discretization data @@ -171,10 +163,19 @@ class LayerPotentialSourceBase(PotentialSource): .. attribute:: complex_dtype .. attribute:: h_max + .. rubric:: Discretization functionality + + See also :class:`meshmode.discretization.Discretization`. + + .. method:: empty + .. method:: nodes + .. method:: num_reference_derivative + .. method:: quad_weights + .. rubric:: Execution - .. method:: weights_and_area_elements .. method:: exec_compute_potential_insn + .. method:: weights_and_area_elements """ @property @@ -197,8 +198,17 @@ class LayerPotentialSourceBase(PotentialSource): def complex_dtype(self): return self.density_discr.complex_dtype + def empty(self, *args, **kwargs): + return self.density_discr.empty(*args, **kwargs) + + def num_reference_derivative(self, *args, **kwargs): + return self.density_discr.num_reference_derivative(*args, **kwargs) + + def quad_weights(self, queue): + return self.density_discr.quad_weights(queue) + @memoize_method - def get_p2p(self, kernels): + def get_p2p(self, kernels, exclude_self=False): # needs to be separate method for caching from pytools import any @@ -209,7 +219,7 @@ class LayerPotentialSourceBase(PotentialSource): from sumpy.p2p import P2P p2p = P2P(self.cl_context, - kernels, exclude_self=False, value_dtypes=value_dtype) + kernels, exclude_self=exclude_self, value_dtypes=value_dtype) return p2p diff --git a/pytential/symbolic/compiler.py b/pytential/symbolic/compiler.py index 68b39d2d..14de3351 100644 --- a/pytential/symbolic/compiler.py +++ b/pytential/symbolic/compiler.py @@ -130,6 +130,8 @@ class PotentialOutput(Record): .. attribute:: kernel_index + .. attribute:: kernel_diagonal_index + .. attribute:: target_name .. attribute:: qbx_forced_limit @@ -139,6 +141,10 @@ class PotentialOutput(Record): center at all) is acceptable. """ + @property + def uses_kernel_diagonal(self): + return self.kernel_diagonal_index is not None + class ComputePotentialInstruction(Instruction): """ @@ -152,6 +158,11 @@ class ComputePotentialInstruction(Instruction): a list of :class:`sumpy.kernel.Kernel` instances, indexed by :attr:`LayerPotentialOutput.kernel_index`. + .. attribute:: kernel_diagonals + + a list of diagonal values to use for the kernel, indexed by + :attr:`LayerPotentialOutput.kernel_diagonal_index` + .. attribute:: kernel_arguments a dictionary mapping arg names to kernel arguments @@ -579,13 +590,21 @@ class OperatorCompiler(IdentityMapper): group = self.group_to_operators[self.op_group_features(expr)] names = [self.get_var_name() for op in group] - kernels = sorted( - set(op.kernel for op in group), - key=lambda kernel: repr(kernel)) + kernels = sorted(set(op.kernel for op in group), key=repr) kernel_to_index = dict( (kernel, i) for i, kernel in enumerate(kernels)) + op_kernel_diagonals = ( + set(op.method_params.get("kernel_diagonal") for op in group) + - set([None])) + + kernel_diagonals = sorted(op_kernel_diagonals, key=repr) + + kernel_diagonal_to_index = dict( + (kernel_diagonal, i) + for i, kernel in enumerate(kernel_diagonals)) + from pytools import single_valued from sumpy.kernel import AxisTargetDerivativeRemover atdr = AxisTargetDerivativeRemover() @@ -603,6 +622,9 @@ class OperatorCompiler(IdentityMapper): PotentialOutput( name=name, kernel_index=kernel_to_index[op.kernel], + kernel_diagonal_index=( + kernel_diagonal_to_index.get( + op.method_params.get("kernel_diagonal"))), target_name=op.target, qbx_forced_limit=op.qbx_forced_limit, ) @@ -613,6 +635,7 @@ class OperatorCompiler(IdentityMapper): ComputePotentialInstruction( outputs=outputs, kernels=tuple(kernels), + kernel_diagonals=tuple(kernel_diagonals), kernel_arguments=kernel_arguments, base_kernel=base_kernel, density=density_var, diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index cca6b919..a3f3cef2 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -365,7 +365,10 @@ class UnregularizedPreprocessor(IdentityMapper): kernel_arguments=dict( (name, self.rec(arg_expr)) for name, arg_expr in expr.kernel_arguments.items() - )) + ), + method_params=dict( + (name, self.rec(arg_expr)) + for name, arg_expr in expr.kernel_arguments.items())) return expr diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index d981fdd2..585f2822 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -672,7 +672,7 @@ class IntG(Expression): def __init__(self, kernel, density, qbx_forced_limit, source=None, target=None, - kernel_arguments=None, + kernel_arguments=None, method_params=None, **kwargs): """*target_derivatives* and later arguments should be considered keyword-only. @@ -680,6 +680,7 @@ class IntG(Expression): :arg kernel: a kernel as accepted by :func:`sumpy.kernel.to_kernel_and_args`, likely a :class:`sumpy.kernel.Kernel`. + :arg qbx_forced_limit: +1 if the output is required to originate from a QBX center on the "+" side of the boundary. -1 for the other side. Evaluation at a target with a value of +/- 1 in *qbx_forced_limit* @@ -696,6 +697,11 @@ class IntG(Expression): ``'avg'`` may be used as a shorthand to evaluate this potential as an average of the ``+1`` and the ``-1`` value. + :arg method_params: Either *None*, or a dictionary containing + method-specific options for evaluation. These include the following: + + * *kernel_diagonal*: The diagonal :math:`K(x,x)` of the kernel + :arg kernel_arguments: A dictionary mapping named :class:`sumpy.kernel.Kernel` arguments (see :meth:`sumpy.kernel.Kernel.get_args` @@ -711,7 +717,10 @@ class IntG(Expression): *kwargs* has the same meaning as *kernel_arguments* can be used as a more user-friendly interface. + """ + if method_params is None: + method_params = {} if kernel_arguments is None: kernel_arguments = {} @@ -771,9 +780,10 @@ class IntG(Expression): self.source = source self.target = target self.kernel_arguments = kernel_arguments + self.method_params = method_params def copy(self, kernel=None, density=None, qbx_forced_limit=_NoArgSentinel, - source=None, target=None, kernel_arguments=None): + source=None, target=None, kernel_arguments=None, method_params=None): kernel = kernel or self.kernel density = density or self.density if qbx_forced_limit is _NoArgSentinel: @@ -781,8 +791,9 @@ class IntG(Expression): source = source or self.source target = target or self.target kernel_arguments = kernel_arguments or self.kernel_arguments + method_params = method_params or self.method_params return type(self)(kernel, density, qbx_forced_limit, source, target, - kernel_arguments) + kernel_arguments, method_params) def __getinitargs__(self): return (self.kernel, self.density, self.qbx_forced_limit, diff --git a/pytential/unregularized.py b/pytential/unregularized.py index 18e8b65c..0717878e 100644 --- a/pytential/unregularized.py +++ b/pytential/unregularized.py @@ -56,13 +56,21 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): """ def __init__(self, density_discr, + area_element_order=None, fmm_order=False, fmm_level_to_order=None, expansion_factory=None, + fmm_backend="sumpy", # begin undocumented arguments # FIXME default debug=False once everything works debug=True): """ + :arg area_element_order: By default, the area elements are numerically + determined by using the same function space as + *density_discr*. Depending on the quadrature rule, this may not be + useful. Passing a non-*None* value of *area_element_order* will + determine the area element using this order discretization. + :arg fmm_order: `False` for direct calculation. """ self.density_discr = density_discr @@ -85,6 +93,8 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): from sumpy.expansion import DefaultExpansionFactory expansion_factory = DefaultExpansionFactory() self.expansion_factory = expansion_factory + self.fmm_backend = fmm_backend + self.area_element_order = area_element_order self.debug = debug @@ -93,27 +103,49 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): import pytential.symbolic.primitives as p from pytential.symbolic.execution import bind with cl.CommandQueue(self.cl_context) as queue: - # quad_stage2_density_discr is not guaranteed to be usable for - # interpolation/differentiation. Use density_discr to find - # area element instead, then upsample that. + # density_discr is not guaranteed to be usable for + # interpolation/differentiation. Use area_element_density_discr to + # find area element instead, then upsample that. - area_element = self.resampler( + area_element = self._area_element_resampler( queue, bind( - self.density_discr, + self._area_element_density_discr, p.area_element(self.ambient_dim, self.dim) )(queue)) - qweight = bind(self.quad_stage2_density_discr, p.QWeight())(queue) + qweight = bind(self.density_discr, p.QWeight())(queue) return (area_element.with_queue(queue)*qweight).with_queue(None) @property - def quad_stage2_density_discr(self): - return self.density_discr + @memoize_method + def _area_element_density_discr(self): + if self.area_element_order is None: + return self.density_discr + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory) - def resampler(self, queue, f): - return f + return Discretization( + self.density_discr.cl_context, self.density_discr.mesh, + InterpolatoryQuadratureSimplexGroupFactory(self.area_element_order), + self.real_dtype) + + @property + @memoize_method + def _area_element_resampler_connection(self): + from meshmode.discretization.connection import make_same_mesh_connection + return make_same_mesh_connection( + self.density_discr, + self._area_element_density_discr) + + def _area_element_resampler(self, queue, f): + if self.area_element_order is None: + return f + + return self._area_element_resampler_connection(queue, f) def with_refinement(self): raise NotImplementedError @@ -132,6 +164,7 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): def exec_compute_potential_insn(self, queue, insn, bound_expr, evaluate): from pytools.obj_array import with_object_array_or_scalar + from functools import partial def evaluate_wrapper(expr): value = evaluate(expr) @@ -175,13 +208,14 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): for o in insn.outputs: target_discr = bound_expr.get_discretization(o.target_name) + is_self = target_discr is self.density_discr if p2p is None: - p2p = self.get_p2p(insn.kernels) + p2p = self.get_p2p(insn.kernels, exclude_self=is_self) evt, output_for_each_kernel = p2p(queue, target_discr.nodes(), - self.quad_stage2_density_discr.nodes(), + self.density_discr.nodes(), [strengths], **kernel_args) result.append((o.name, output_for_each_kernel[o.kernel_index])) diff --git a/test/test_unregularized.py b/test/test_unregularized.py index c104acd4..c7a0efb5 100644 --- a/test/test_unregularized.py +++ b/test/test_unregularized.py @@ -46,8 +46,8 @@ circle = partial(ellipse, 1) # {{{ ones kernel test -def test_ones_kernel(ctx_getter): - cl_ctx = ctx_getter() +def test_ones_kernel(ctx_factory): + cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) nelements = 10 @@ -89,20 +89,77 @@ def test_ones_kernel(ctx_getter): # }}} -# {{{ on surface direct +# {{{ on surface tests -# }}} +def run_on_surface_eigval_test(cl_ctx, fmm_order, fmm_backend): + nelements = 100 + target_order = 16 + quad_order = 0 + + aspect = 2 + mode_nr = 4 + + queue = cl.CommandQueue(cl_ctx) + + mesh = make_curve_mesh( + partial(ellipse, aspect), + np.linspace(0, 1, nelements + 1), + target_order) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + density_discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(quad_order)) + + nodes = density_discr.nodes().with_queue(queue) + + from pytential.unregularized import UnregularizedLayerPotentialSource + lpot_source = UnregularizedLayerPotentialSource( + density_discr, fmm_order=fmm_order, + fmm_backend=fmm_backend, + area_element_order=target_order) + + from sumpy.kernel import LaplaceKernel + knl = LaplaceKernel(2) + + sym_kernel_diagonal = (-1 / (2 * np.pi)) * sym.mean_curvature(2) + + d_op = sym.D(knl, sym.var("sigma"), + method_params=dict(kernel_diagonal=sym_kernel_diagonal), + qbx_forced_limit=None) + + angle = cl.clmath.atan2(nodes[1], nodes[0]) + sigma = cl.clmath.cos(mode_nr * angle) + + ellipse_fraction = ((1-aspect)/(1+aspect))**mode_nr + + bound_d_op = bind(lpot_source, d_op) + + d_sigma = bound_d_op(queue=queue, sigma=sigma) + + d_eigval = -(-1)**mode_nr * 1/2*ellipse_fraction + d_sigma_ref = d_eigval*sigma + d_ref_norm = norm(density_discr, queue, d_sigma_ref) + + print(norm(density_discr, queue, (d_sigma - d_sigma_ref) / d_ref_norm)) -# {{{ on surface fmm +@pytest.mark.parametrize(("fmm_order", "fmm_backend"), [ + (False, None)]) +def test_on_surface_eigvals(ctx_factory, fmm_order, fmm_backend): + cl_ctx = ctx_factory() + run_on_surface_eigval_test(cl_ctx, fmm_order, fmm_backend) # }}} # {{{ off surface fmm versus direct -def test_off_surface_fmm_vs_direct(ctx_getter): - cl_ctx = ctx_getter() +def test_off_surface_fmm_vs_direct(ctx_factory): + cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) nelements = 300 -- GitLab