From 393b876d2428fbc65e0099502e7718472d6b4548 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 4 Apr 2017 00:58:51 -0500 Subject: [PATCH 1/2] Add point potential source, switch IE test to use it --- pytential/qbx/__init__.py | 25 ++---- pytential/source.py | 132 +++++++++++++++++++++++++++++++ pytential/symbolic/compiler.py | 12 +-- pytential/symbolic/execution.py | 4 +- pytential/symbolic/primitives.py | 28 ++++++- test/test_layer_pot.py | 55 ++++++------- 6 files changed, 197 insertions(+), 59 deletions(-) create mode 100644 pytential/source.py diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 71a72ce6..a18196f4 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -31,6 +31,7 @@ from pytools import memoize_method from meshmode.discretization import Discretization from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory from pytential.qbx.target_assoc import QBXTargetAssociationFailedException +from pytential.source import PotentialSource import pyopencl as cl @@ -84,18 +85,8 @@ class _JumpTermArgumentProvider(object): # }}} -class LayerPotentialSource(object): - """ - .. method:: preprocess_optemplate(name, expr) - - .. method:: op_group_features(expr) - - Return a characteristic tuple by which operators that can be - executed together can be grouped. - - *expr* is a subclass of - :class:`pytential.symbolic.primitives.IntG`. - """ +class LayerPotentialSource(PotentialSource): + pass def get_local_expansion_class(base_kernel): @@ -430,7 +421,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): return result - def exec_layer_potential_insn(self, queue, insn, bound_expr, evaluate): + 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 oversample = partial(self.resampler, queue) @@ -446,9 +437,9 @@ class QBXLayerPotentialSource(LayerPotentialSource): return with_object_array_or_scalar(oversample, value) if self.fmm_level_to_order is False: - func = self.exec_layer_potential_insn_direct + func = self.exec_compute_potential_insn_direct else: - func = self.exec_layer_potential_insn_fmm + func = self.exec_compute_potential_insn_fmm return func(queue, insn, bound_expr, evaluate_wrapper) @@ -496,7 +487,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_mpole_factory, fmm_local_factory, qbx_local_factory, out_kernels) - def exec_layer_potential_insn_fmm(self, queue, insn, bound_expr, evaluate): + def exec_compute_potential_insn_fmm(self, queue, insn, bound_expr, evaluate): # {{{ build list of unique target discretizations used # map (name, qbx_side) to number in list @@ -706,7 +697,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): *count = item; """) - def exec_layer_potential_insn_direct(self, queue, insn, bound_expr, evaluate): + def exec_compute_potential_insn_direct(self, queue, insn, bound_expr, evaluate): lpot_applier = self.get_lpot_applier(insn.kernels) p2p = None lpot_applier_on_tgt_subset = None diff --git a/pytential/source.py b/pytential/source.py new file mode 100644 index 00000000..87565060 --- /dev/null +++ b/pytential/source.py @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2017 Andreas Kloeckner" + +__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 # noqa: F401 +import pyopencl as cl # noqa: F401 +import six +from pytools import memoize_method + + +class PotentialSource(object): + """ + .. method:: preprocess_optemplate(name, expr) + + .. method:: op_group_features(expr) + + Return a characteristic tuple by which operators that can be + executed together can be grouped. + + *expr* is a subclass of + :class:`pytential.symbolic.primitives.IntG`. + """ + + +class PointPotentialSource(PotentialSource): + """ + ... attributes:: points + + An :class:`pyopencl.array.Array` of shape ``[ambient_dim, npoints]``. + + """ + + def __init__(self, cl_context, points): + self.cl_context = cl_context + self.points = points + + @property + def real_dtype(self): + return self.points.dtype + + @property + def complex_dtype(self): + return { + np.float32: np.complex64, + np.float64: np.complex128 + }[self.real_dtype.type] + + @property + def ambient_dim(self): + return self.points.shape[0] + + def op_group_features(self, expr): + from sumpy.kernel import AxisTargetDerivativeRemover + result = ( + expr.source, expr.density, + AxisTargetDerivativeRemover()(expr.kernel), + ) + + return result + + @memoize_method + def get_p2p(self, kernels): + # needs to be separate method for caching + + from pytools import any + if any(knl.is_complex_valued for knl in kernels): + value_dtype = self.complex_dtype + else: + value_dtype = self.real_dtype + + from sumpy.p2p import P2P + p2p = P2P(self.cl_context, + kernels, exclude_self=False, value_dtypes=value_dtype) + + return p2p + + def exec_compute_potential_insn(self, queue, insn, bound_expr, evaluate): + p2p = None + + kernel_args = {} + for arg_name, arg_expr in six.iteritems(insn.kernel_arguments): + kernel_args[arg_name] = evaluate(arg_expr) + + strengths = (evaluate(insn.density).with_queue(queue) + * self.weights_and_area_elements()) + + # FIXME: Do this all at once + result = [] + for o in insn.outputs: + target_discr = bound_expr.get_discretization(o.target_name) + + # no on-disk kernel caching + if p2p is None: + p2p = self.get_p2p(insn.kernels) + + evt, output_for_each_kernel = p2p(queue, + target_discr.nodes(), self.points, + [strengths], **kernel_args) + + result.append((o.name, output_for_each_kernel[o.kernel_index])) + + return result, [] + + @memoize_method + def weights_and_area_elements(self): + with cl.CommandQueue(self.cl_context) as queue: + result = cl.array.empty(queue, self.points.shape[-1], + dtype=self.real_dtype) + result.fill(1) + + return result.with_queue(None) diff --git a/pytential/symbolic/compiler.py b/pytential/symbolic/compiler.py index 31379ce3..fc97155f 100644 --- a/pytential/symbolic/compiler.py +++ b/pytential/symbolic/compiler.py @@ -122,7 +122,7 @@ class Assign(Instruction): # {{{ layer pot instruction -class LayerPotentialOutput(Record): +class PotentialOutput(Record): """ .. attribute:: name @@ -140,11 +140,11 @@ class LayerPotentialOutput(Record): """ -class LayerPotentialInstruction(Instruction): +class ComputePotentialInstruction(Instruction): """ .. attribute:: outputs - A list of :class:`LayerPotentialOutput` instances + A list of :class:`PotentialOutput` instances The entries in the list correspond to :attr:`names`. .. attribute:: kernels @@ -225,7 +225,7 @@ class LayerPotentialInstruction(Instruction): def get_exec_function(self, exec_mapper): source = exec_mapper.bound_expr.places[self.source] - return source.exec_layer_potential_insn + return source.exec_compute_potential_insn def __hash__(self): return id(self) @@ -600,7 +600,7 @@ class OperatorCompiler(IdentityMapper): for arg_name, arg_val in six.iteritems(expr.kernel_arguments)) outputs = [ - LayerPotentialOutput( + PotentialOutput( name=name, kernel_index=kernel_to_index[op.kernel], target_name=op.target, @@ -610,7 +610,7 @@ class OperatorCompiler(IdentityMapper): ] self.code.append( - LayerPotentialInstruction( + ComputePotentialInstruction( outputs=outputs, kernels=tuple(kernels), kernel_arguments=kernel_arguments, diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index e2847243..a229b48e 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -326,8 +326,8 @@ def prepare_places(places): def cast_to_place(discr): from pytential.target import TargetBase - if not isinstance(discr, (Discretization, TargetBase, - LayerPotentialSource)): + from pytential.source import PotentialSource + if not isinstance(discr, (Discretization, TargetBase, PotentialSource)): raise TypeError("must pass discretizations, " "layer potential sources or targets as 'places'") return discr diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 0a2d28ea..27c2e24c 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -27,14 +27,17 @@ from six.moves import intern from warnings import warn import numpy as np -from pymbolic.primitives import ( # noqa +from pymbolic.primitives import ( # noqa: F401, N813 Expression as ExpressionBase, Variable as var, cse_scope as cse_scope_base, make_common_subexpression as cse) from pymbolic.geometric_algebra import MultiVector, componentwise -from pymbolic.geometric_algebra.primitives import ( # noqa +from pymbolic.geometric_algebra.primitives import ( # noqa: F401 NablaComponent, DerivativeSource, Derivative as DerivativeBase) -from pymbolic.primitives import make_sym_vector # noqa +from pymbolic.primitives import make_sym_vector # noqa: F401 +from pytools.obj_array import make_obj_array # noqa: F401 + +from functools import partial __doc__ = """ @@ -203,7 +206,6 @@ def nodes(ambient_dim, where=None): locations. """ - from pytools.obj_array import make_obj_array return MultiVector( make_obj_array([ NodeCoordinateComponent(i, where) @@ -415,6 +417,24 @@ class Derivative(DerivativeBase): return DerivativeBinder()(expr) +def dd_axis(axis, ambient_dim, operand): + d = Derivative() + + unit_vector = np.zeros(ambient_dim) + unit_vector[axis] = 1 + + unit_mvector = MultiVector(unit_vector) + + return d.resolve( + (unit_mvector.scalar_product(d.dnabla(ambient_dim))) + * d(operand)) + + +d_dx = partial(dd_axis, 0) +d_dy = partial(dd_axis, 1) +d_dz = partial(dd_axis, 2) + + # {{{ potentials def hashable_kernel_args(kernel_arguments): diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 1c3f2598..87c11992 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -322,7 +322,7 @@ def run_int_eq_test( DirichletOperator, NeumannOperator) - from sumpy.kernel import LaplaceKernel, HelmholtzKernel, AxisTargetDerivative + from sumpy.kernel import LaplaceKernel, HelmholtzKernel if k: knl = HelmholtzKernel(2) knl_kwargs = {"k": sym.var("k")} @@ -372,6 +372,8 @@ def run_int_eq_test( source_charges = source_charges.astype(dtype) assert np.sum(source_charges) < 1e-15 + source_charges_dev = cl.array.to_device(queue, source_charges) + # }}} if 0: @@ -385,32 +387,28 @@ def run_int_eq_test( # {{{ establish BCs - from sumpy.p2p import P2P - pot_p2p = P2P(cl_ctx, - [knl], exclude_self=False, value_dtypes=dtype) - - evt, (test_direct,) = pot_p2p( - queue, test_targets, point_sources, [source_charges], - out_host=False, **concrete_knl_kwargs) + from pytential.source import PointPotentialSource + from pytential.target import PointsTarget - nodes = density_discr.nodes() + point_source = PointPotentialSource(cl_ctx, point_sources) - evt, (src_pot,) = pot_p2p( - queue, nodes, point_sources, [source_charges], - **concrete_knl_kwargs) + pot_src = sym.IntG( + # FIXME: qbx_forced_limit--really? + knl, sym.var("charges"), qbx_forced_limit=None, k=sym.var("k")) - grad_p2p = P2P(cl_ctx, - [AxisTargetDerivative(0, knl), AxisTargetDerivative(1, knl)], - exclude_self=False, value_dtypes=dtype) - evt, (src_grad0, src_grad1) = grad_p2p( - queue, nodes, point_sources, [source_charges], - **concrete_knl_kwargs) + test_direct = bind((point_source, PointsTarget(test_targets)), pot_src)( + queue, charges=source_charges_dev, **concrete_knl_kwargs) if bc_type == "dirichlet": - bc = src_pot + bc = bind((point_source, density_discr), pot_src)( + queue, charges=source_charges_dev, **concrete_knl_kwargs) + elif bc_type == "neumann": - normal = bind(density_discr, sym.normal(2))(queue).as_vector(np.object) - bc = (src_grad0*normal[0] + src_grad1*normal[1]) + bc = bind( + (point_source, density_discr), + sym.normal_derivative( + qbx.ambient_dim, pot_src, where=sym.DEFAULT_TARGET) + )(queue, charges=source_charges_dev, **concrete_knl_kwargs) # }}} @@ -450,8 +448,6 @@ def run_int_eq_test( # {{{ error check - from pytential.target import PointsTarget - bound_tgt_op = bind((qbx, PointsTarget(test_targets)), op.representation(sym.var("u"))) @@ -493,19 +489,18 @@ def run_int_eq_test( tang_deriv_from_src = bound_t_deriv_op( queue, u=u, **concrete_knl_kwargs).as_scalar().get() - tangent = bind( - density_discr, - sym.pseudoscalar(2)/sym.area_element(2))( - queue, **concrete_knl_kwargs).as_vector(np.object) - - tang_deriv_ref = (src_grad0 * tangent[0] + src_grad1 * tangent[1]).get() + tang_deriv_ref = (bind( + (point_source, density_discr), + sym.tangential_derivative(2, pot_src) + )(queue, charges=source_charges_dev, **concrete_knl_kwargs) + .as_scalar().get()) if 0: pt.plot(tang_deriv_ref.real) pt.plot(tang_deriv_from_src.real) pt.show() - td_err = tang_deriv_from_src - tang_deriv_ref + td_err = (tang_deriv_from_src - tang_deriv_ref) rel_td_err_inf = la.norm(td_err, np.inf)/la.norm(tang_deriv_ref, np.inf) -- GitLab From 34b403704be8e7092fc3e4af8d9b33713d6e41e4 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 4 Apr 2017 01:53:52 -0500 Subject: [PATCH 2/2] More test_layer_pot fixes for point pot source --- test/test_layer_pot.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 87c11992..d89615ec 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -394,7 +394,7 @@ def run_int_eq_test( pot_src = sym.IntG( # FIXME: qbx_forced_limit--really? - knl, sym.var("charges"), qbx_forced_limit=None, k=sym.var("k")) + knl, sym.var("charges"), qbx_forced_limit=None, **knl_kwargs) test_direct = bind((point_source, PointsTarget(test_targets)), pot_src)( queue, charges=source_charges_dev, **concrete_knl_kwargs) @@ -518,9 +518,8 @@ def run_int_eq_test( #pt.plot(u) #pt.show() - evt, (fld_from_src,) = pot_p2p( - queue, fplot.points, point_sources, [source_charges], - **knl_kwargs) + fld_from_src = bind((point_source, PointsTarget(fplot.points)), + pot_src)(queue, charges=source_charges_dev, **concrete_knl_kwargs) fld_from_bdry = bind( (qbx, PointsTarget(fplot.points)), op.representation(sym.var("u")) -- GitLab