diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 2d85ee1036117363cccffdf5c2eee81da57c4fcc..bc29d7faec9979b7c9555aab2fd90e68e46947b4 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -618,31 +618,18 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): # {{{ internal functionality for execution def exec_compute_potential_insn(self, queue, insn, bound_expr, evaluate): - from pytools.obj_array import with_object_array_or_scalar - - def oversample_nonscalars(vec): - from numbers import Number - if isinstance(vec, Number): - return vec - else: - return self.resampler(queue, vec) - if not self._refined_for_global_qbx: from warnings import warn warn( "Executing global QBX without refinement. " "This is unlikely to work.") - def evaluate_wrapper(expr): - value = evaluate(expr) - return with_object_array_or_scalar(oversample_nonscalars, value) - if self.fmm_level_to_order is False: func = self.exec_compute_potential_insn_direct else: func = self.exec_compute_potential_insn_fmm - return func(queue, insn, bound_expr, evaluate_wrapper) + return func(queue, insn, bound_expr, evaluate) @property @memoize_method diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 165ab796633e9deaf352e29521513a54e7fa7366..9a9bcb05dd2cf18ce70ce9f11bfb448be3ef2d69 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -39,9 +39,7 @@ import pyopencl.clmath # noqa from loopy.version import MOST_RECENT_LANGUAGE_VERSION from pytools import memoize_in -from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET -from pytential.symbolic.primitives import ( - QBXSourceStage1, QBXSourceStage2, QBXSourceQuadStage2) +from pytential import sym # FIXME caches: fix up queues @@ -175,6 +173,41 @@ class EvaluationMapper(EvaluationMapperBase): source = self.bound_expr.places[expr.source] return source.map_quad_kernel_op(expr, self.bound_expr, self.rec) + def map_interpolation(self, expr): + if not isinstance(expr.target, sym.QBXSourceQuadStage2): + raise RuntimeError("can only interpolate to quad_stage2 mesh") + + source = self.bound_expr.places[expr.source] + operand = self.rec(expr.operand) + + from pytential.source import LayerPotentialSourceBase + if isinstance(source, LayerPotentialSourceBase): + sym_source = expr.source + if not isinstance(sym_source, sym._QBXSource): + sym_source = sym.QBXSourceStage1(sym_source) + + if isinstance(sym_source, sym.QBXSourceStage1): + resampler = source.resampler + elif isinstance(sym_source, sym.QBXSourceStage2): + resampler = source.refined_interp_to_ovsmp_quad_connection + elif isinstance(sym_source, sym.QBXSourceQuadStage2): + resampler = lambda x, y: y # noqa + else: + from pytential.symbolic.mappers import stringify_where + raise ValueError( + "unknown `where` identifier in " + "interpolation source: {}".format( + stringify_where(sym_source))) + else: + raise TypeError("source must be a `LayerPotentialSourceBase`") + + if isinstance(operand, cl.array.Array): + return resampler(self.queue, operand) + elif isinstance(operand, (int, float, complex, np.number)): + return operand + else: + raise TypeError("cannot interpolate `{}`".format(type(operand))) + # }}} def exec_assign(self, queue, insn, bound_expr, evaluate): @@ -375,7 +408,7 @@ class GeometryCollection(object): from pytential.source import LayerPotentialSourceBase, PotentialSource if auto_where is None: - source_where, target_where = DEFAULT_SOURCE, DEFAULT_TARGET + source_where, target_where = sym.DEFAULT_SOURCE, sym.DEFAULT_TARGET else: # NOTE: keeping this here to make sure auto_where unpacks into # just the two elements @@ -389,7 +422,7 @@ class GeometryCollection(object): if isinstance(places, LayerPotentialSourceBase): self.places[source_where] = places self.places[target_where] = \ - self._get_lpot_discretization(target_where, places) + self._get_lpot_discretization(places, target_where) elif isinstance(places, (Discretization, TargetBase)): self.places[target_where] = places elif isinstance(places, tuple): @@ -406,23 +439,20 @@ class GeometryCollection(object): self.caches = {} - def _get_lpot_discretization(self, where, lpot): - from pytential.source import LayerPotentialSourceBase - if not isinstance(lpot, LayerPotentialSourceBase): - return lpot + def _get_lpot_discretization(self, lpot, where): + if not isinstance(where, sym._QBXSource): + where = sym.QBXSourceStage1(where) - from pytential.symbolic.primitives import _QBXSource - if not isinstance(where, _QBXSource): - where = QBXSourceStage1(where) - - if isinstance(where, QBXSourceStage1): + if isinstance(where, sym.QBXSourceStage1): return lpot.density_discr - if isinstance(where, QBXSourceStage2): + if isinstance(where, sym.QBXSourceStage2): return lpot.stage2_density_discr - if isinstance(where, QBXSourceQuadStage2): + if isinstance(where, sym.QBXSourceQuadStage2): return lpot.quad_stage2_density_discr - raise ValueError('unknown `where` identifier: {}'.format(type(where))) + from pytential.symbolic.mappers import stringify_where + raise ValueError('unknown `where` identifier: {}'.format( + stringify_where(where))) def get_discretization(self, where): """ @@ -436,14 +466,18 @@ class GeometryCollection(object): """ if where in self.places: - lpot = self.places[where] + discr = self.places[where] else: - lpot = self.places.get(getattr(where, 'where', None), None) + discr = self.places.get(getattr(where, 'where', None), None) - if lpot is None: + if discr is None: raise KeyError('`where` not in the collection: {}'.format(where)) - return self._get_lpot_discretization(where, lpot) + from pytential.source import LayerPotentialSourceBase + if isinstance(discr, LayerPotentialSourceBase): + return self._get_lpot_discretization(discr, where) + else: + return discr def __getitem__(self, where): return self.places[where] @@ -492,7 +526,7 @@ class BoundExpression(object): nresults = 1 domains = _prepare_domains(nresults, self.places, domains, - DEFAULT_TARGET) + sym.DEFAULT_TARGET) total_dofs = 0 starts_and_ends = [] @@ -530,7 +564,9 @@ def bind(places, expr, auto_where=None): if not isinstance(places, GeometryCollection): places = GeometryCollection(places, auto_where=auto_where) + from pytential.symbolic.mappers import QBXInterpolationPreprocessor expr = _prepare_expr(places, expr) + expr = QBXInterpolationPreprocessor(places)(expr) return BoundExpression(places, expr) diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 00c43689e45ffbe435a87f9b8e6fd3c497d87303..415453537e877bb40d7dbe35ade475e0337e51ef 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -103,6 +103,9 @@ class IdentityMapper(IdentityMapperBase): for name, arg_expr in expr.kernel_arguments.items() )) + def map_interpolation(self, expr): + return type(expr)(expr.source, expr.target, self.rec(expr.operand)) + class CombineMapper(CombineMapperBase): def map_node_sum(self, expr): @@ -113,6 +116,7 @@ class CombineMapper(CombineMapperBase): map_elementwise_sum = map_node_sum map_elementwise_min = map_node_sum map_elementwise_max = map_node_sum + map_interpolation = map_node_sum def map_int_g(self, expr): return self.combine( @@ -280,6 +284,17 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper): for name, name_expr in six.iteritems(expr.extra_vars)]), where) + def map_interpolation(self, expr): + source = expr.source + target = expr.target + + if source is None: + source = self.default_source + if target is None: + target = self.default_where + + return type(expr)(source, target, expr.operand) + def operand_rec(self, expr): return self.rec(expr) @@ -393,6 +408,35 @@ class UnregularizedPreprocessor(IdentityMapper): # {{{ QBX preprocessor +class QBXInterpolationPreprocessor(IdentityMapper): + def __init__(self, places): + self.places = places + + def map_int_g(self, expr): + if isinstance(expr.source, prim.QBXSourceQuadStage2): + return expr + + from pytential.source import LayerPotentialSourceBase + source = self.places[expr.source] + + if isinstance(source, LayerPotentialSourceBase): + stage2_source = prim.QBXSourceQuadStage2(prim.DEFAULT_SOURCE) + + density = prim.Interpolation( + expr.source, stage2_source, self.rec(expr.density)) + kernel_arguments = dict( + (name, prim.Interpolation( + expr.source, stage2_source, self.rec(arg_expr))) + for name, arg_expr in expr.kernel_arguments.items()) + + expr = expr.copy( + kernel=expr.kernel, + density=density, + kernel_arguments=kernel_arguments) + + return expr + + class QBXPreprocessor(IdentityMapper): def __init__(self, source_name, places): self.source_name = source_name @@ -554,6 +598,12 @@ class StringifyMapper(BaseStringifyMapper): expr.kernel, self.rec(expr.density, PREC_PRODUCT)) + def map_interpolation(self, expr, enclosing_prec): + return "Interp[%s->%s](%s)" % ( + stringify_where(expr.source), + stringify_where(expr.target), + self.rec(expr.operand, PREC_PRODUCT)) + # }}} diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index f7332a450f593ff094640d87af9556e68754e422..ec129369422f1a04d3f16144df47240b739269b7 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -141,6 +141,11 @@ Elementary numerics .. autofunction:: mean .. autoclass:: IterativeInverse +Operators +^^^^^^^^^ + +.. autoclass:: Interpolation + Geometric Calculus (based on Geometric/Clifford Algebra) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -1017,6 +1022,31 @@ class _NoArgSentinel(object): pass +class Interpolation(Expression): + """Interpolate quantity from *source* to *target* discretization.""" + + init_arg_names = ("source", "target", "operand") + + def __new__(cls, source, target, operand): + if isinstance(operand, np.ndarray): + def make_op(operand_i): + return cls(source, target, operand_i) + + return componentwise(make_op, operand) + else: + return Expression.__new__(cls) + + def __init__(self, source, target, operand): + self.source = source + self.target = target + self.operand = operand + + def __getinitargs__(self): + return (self.source, self.target, self.operand) + + mapper_method = intern("map_interpolation") + + class IntG(Expression): r""" .. math:: diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py index 1e048b7f6495e664767850c18b6e9bf7f6c2d404..b65b5263f20173c5911c1077c2ae46e6d8bd01cc 100644 --- a/test/test_scalar_int_eq.py +++ b/test/test_scalar_int_eq.py @@ -630,7 +630,6 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize): # {{{ solve bound_op = bind(qbx, op_u) - rhs = bind(density_discr, op.prepare_rhs(sym.var("bc")))(queue, bc=bc) try: diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 22a2a2fd23dc792951a10364c291406b09af11e7..d709f7f0e57cebe40bdf38be90e48e39c8403e5e 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -26,6 +26,7 @@ import pytest import numpy as np import numpy.linalg as la + import pyopencl as cl import pyopencl.array # noqa import pyopencl.clmath # noqa @@ -176,15 +177,13 @@ def test_tangential_onb(ctx_factory): # {{{ test_expr_pickling def test_expr_pickling(): - from sumpy.kernel import LaplaceKernel, AxisTargetDerivative import pickle + from sumpy.kernel import LaplaceKernel, AxisTargetDerivative ops_for_testing = [ sym.d_dx( 2, - sym.D( - LaplaceKernel(2), sym.var("sigma"), qbx_forced_limit=-2 - ) + sym.D(LaplaceKernel(2), sym.var("sigma"), qbx_forced_limit=-2) ), sym.D( AxisTargetDerivative(0, LaplaceKernel(2)), @@ -202,6 +201,53 @@ def test_expr_pickling(): # }}} +@pytest.mark.parametrize("source", [ + sym.DEFAULT_SOURCE, + sym.QBXSourceStage1(sym.DEFAULT_SOURCE), + sym.QBXSourceStage2(sym.DEFAULT_SOURCE), + sym.QBXSourceQuadStage2(sym.DEFAULT_SOURCE) + ]) +def test_interpolation(ctx_factory, source): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + nelements = 32 + target_order = 7 + qbx_order = 4 + + mesh = make_curve_mesh(starfish, + np.linspace(0.0, 1.0, nelements + 1), + target_order) + discr = Discretization(ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + from pytential.qbx import QBXLayerPotentialSource + qbx, _ = QBXLayerPotentialSource(discr, + fine_order=4 * target_order, + qbx_order=qbx_order, + fmm_order=False).with_refinement() + + target = sym.QBXSourceQuadStage2(sym.DEFAULT_SOURCE) + sigma_sym = sym.var("sigma") + op_sym = sym.sin(sym.Interpolation(source, target, sigma_sym)) + bound_op = bind(qbx, op_sym, auto_where=(source, sym.DEFAULT_TARGET)) + + quad2_nodes = qbx.quad_stage2_density_discr.nodes().get(queue) + if isinstance(source, sym.QBXSourceStage2): + nodes = qbx.stage2_density_discr.nodes().get(queue) + elif isinstance(source, sym.QBXSourceQuadStage2): + nodes = quad2_nodes + else: + nodes = qbx.density_discr.nodes().get(queue) + + sigma_dev = cl.array.to_device(queue, la.norm(nodes, axis=0)) + sigma_quad2 = np.sin(la.norm(quad2_nodes, axis=0)) + sigma_quad2_interp = bound_op(queue, sigma=sigma_dev).get(queue) + + error = la.norm(sigma_quad2_interp - sigma_quad2) / la.norm(sigma_quad2) + assert error < 1.0e-10 + + # You can test individual routines by typing # $ python test_symbolic.py 'test_routine()'