From bfdee409657044b23970d6c834a2fb28c9aa463b Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Wed, 7 Nov 2018 22:03:23 -0600 Subject: [PATCH 01/11] start adding an interpolation preprocessor --- pytential/symbolic/execution.py | 21 +++++++++++++++++++++ pytential/symbolic/mappers.py | 9 +++++++++ pytential/symbolic/primitives.py | 15 +++++++++++++++ 3 files changed, 45 insertions(+) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 165ab796..f949fb66 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -175,6 +175,27 @@ 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): + assert isinstance(expr.target, sym.QBXSourceQuadStage2) + + where = expr.source + if isinstance(where, sym.DEFAULT_SOURCE): + where = sym.QBXSourceStage1(where) + + source = self.bound_expr.places[expr.source] + density = self.rec(expr.density) + + if isinstance(where, sym.QBXSourceStage1): + density = source.resampler(density) + elif isinstance(where, sym.QBXSourceStage2): + density = source.refined_interp_to_ovsmp_quad_connection(density) + elif isinstance(where, sym.QBXSourceQuadStage2): + pass + else: + raise ValueError('unknown `where` identifier: {}'.format(type(where))) + + return density + # }}} def exec_assign(self, queue, insn, bound_expr, evaluate): diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 00c43689..0c930f6f 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -393,6 +393,15 @@ class UnregularizedPreprocessor(IdentityMapper): # {{{ QBX preprocessor +class InterpolationPreprocessor(IdentityMapper): + def mat_int_g(self, expr): + if isinstance(expr.source, prim.QBXSourceQuadStage2): + return expr + + density = prim.Interpolation(self.rec(expr.density), expr.source) + return expr.copy(density=density) + + class QBXPreprocessor(IdentityMapper): def __init__(self, source_name, places): self.source_name = source_name diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index f7332a45..833eff0d 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -1017,6 +1017,21 @@ class _NoArgSentinel(object): pass +class Interpolation(Expression): + """Interpolate density from *source* to *target* discretization.""" + + init_arg_names = ("density", "source", "target") + + def __init__(self, density, source, target=None): + self.density = density + self.source = source + self.target = target + if self.target is None: + self.target = QBXSourceQuadStage2(DEFAULT_SOURCE) + + mapper_method = intern("map_interpolation") + + class IntG(Expression): r""" .. math:: -- GitLab From ae22246a801b6b905b294c403eccd0ae16da56b2 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Thu, 8 Nov 2018 20:47:22 -0600 Subject: [PATCH 02/11] hook up mappers to the interpolation primitive --- pytential/symbolic/execution.py | 40 +++++++++++------------ pytential/symbolic/mappers.py | 15 +++++++++ pytential/symbolic/primitives.py | 3 ++ test/test_symbolic.py | 54 +++++++++++++++++++++++++++++--- 4 files changed, 88 insertions(+), 24 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index f949fb66..65c1dcc7 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 @@ -179,22 +177,23 @@ class EvaluationMapper(EvaluationMapperBase): assert isinstance(expr.target, sym.QBXSourceQuadStage2) where = expr.source - if isinstance(where, sym.DEFAULT_SOURCE): + if not isinstance(where, sym._QBXSource): where = sym.QBXSourceStage1(where) source = self.bound_expr.places[expr.source] - density = self.rec(expr.density) - if isinstance(where, sym.QBXSourceStage1): - density = source.resampler(density) + resampler = source.resampler elif isinstance(where, sym.QBXSourceStage2): - density = source.refined_interp_to_ovsmp_quad_connection(density) + resampler = source.refined_interp_to_ovsmp_quad_connection elif isinstance(where, sym.QBXSourceQuadStage2): - pass + resampler = lambda x, y: y else: - raise ValueError('unknown `where` identifier: {}'.format(type(where))) + from pytential.symbolic.mappers import stringify_where + raise ValueError('unknown `where` identifier: {}'.format( + stringify_where(where))) - return density + density = self.rec(expr.density) + return resampler(self.queue, density) # }}} @@ -396,7 +395,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 @@ -432,18 +431,19 @@ class GeometryCollection(object): if not isinstance(lpot, LayerPotentialSourceBase): return lpot - from pytential.symbolic.primitives import _QBXSource - if not isinstance(where, _QBXSource): - where = QBXSourceStage1(where) + if not isinstance(where, sym._QBXSource): + where = sym.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( + sym.stringify_where(where))) def get_discretization(self, where): """ @@ -513,7 +513,7 @@ class BoundExpression(object): nresults = 1 domains = _prepare_domains(nresults, self.places, domains, - DEFAULT_TARGET) + sym.DEFAULT_TARGET) total_dofs = 0 starts_and_ends = [] diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 0c930f6f..79b641ba 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -81,6 +81,7 @@ class IdentityMapper(IdentityMapperBase): map_node_coordinate_component = map_ones map_parametrization_gradient = map_ones map_parametrization_derivative = map_ones + map_interpolation = map_ones # }}} @@ -108,6 +109,9 @@ class CombineMapper(CombineMapperBase): def map_node_sum(self, expr): return self.rec(expr.operand) + def map_interpolation(self, expr): + return self.rec(expr.density) + map_node_max = map_node_sum map_num_reference_derivative = map_node_sum map_elementwise_sum = map_node_sum @@ -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)(expr.density, source, target) + def operand_rec(self, expr): return self.rec(expr) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 833eff0d..a19d97c8 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -1029,6 +1029,9 @@ class Interpolation(Expression): if self.target is None: self.target = QBXSourceQuadStage2(DEFAULT_SOURCE) + def __getinitargs__(self): + return (self.density, self.source, self.target) + mapper_method = intern("map_interpolation") diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 22a2a2fd..33efe9f9 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("where", [ + sym.DEFAULT_SOURCE, + sym.QBXSourceStage1(sym.DEFAULT_SOURCE), + sym.QBXSourceStage2(sym.DEFAULT_SOURCE), + sym.QBXSourceQuadStage2(sym.DEFAULT_SOURCE) + ]) +def test_interpolation(ctx_factory, where): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + ndim = 2 + nelements = 32 + target_order = 5 + 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() + + sigma_sym = sym.var("sigma") + op_sym = sym.sin(sym.Interpolation(sigma_sym, where)) + bound_op = bind(qbx, op_sym, auto_where=(where, sym.DEFAULT_TARGET)) + + quad2_nodes = qbx.quad_stage2_density_discr.nodes().get(queue) + if isinstance(where, sym.QBXSourceStage2): + nodes = qbx.stage2_density_discr.nodes().get(queue) + elif isinstance(where, 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) + print(error) + + # You can test individual routines by typing # $ python test_symbolic.py 'test_routine()' -- GitLab From 33896b0bf5b5c62f413090376cbbc5361bc2aa65 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sat, 10 Nov 2018 20:51:52 -0600 Subject: [PATCH 03/11] hook up interpolation into bind --- pytential/qbx/__init__.py | 15 +-------------- pytential/symbolic/execution.py | 19 ++++++++++++++----- pytential/symbolic/mappers.py | 27 ++++++++++++++++++++------- pytential/symbolic/primitives.py | 19 ++++++++++++++----- test/test_scalar_int_eq.py | 1 - test/test_symbolic.py | 3 +-- 6 files changed, 50 insertions(+), 34 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 2d85ee10..bc29d7fa 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 65c1dcc7..0d60ac29 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -176,24 +176,32 @@ class EvaluationMapper(EvaluationMapperBase): def map_interpolation(self, expr): assert isinstance(expr.target, sym.QBXSourceQuadStage2) + source = self.bound_expr.places[expr.source] + operand = self.rec(expr.operand) + + from pytential.source import LayerPotentialSourceBase + if not isinstance(source, LayerPotentialSourceBase): + return operand + + if not isinstance(operand, cl.array.Array): + return operand + where = expr.source if not isinstance(where, sym._QBXSource): where = sym.QBXSourceStage1(where) - source = self.bound_expr.places[expr.source] if isinstance(where, sym.QBXSourceStage1): resampler = source.resampler elif isinstance(where, sym.QBXSourceStage2): resampler = source.refined_interp_to_ovsmp_quad_connection elif isinstance(where, sym.QBXSourceQuadStage2): - resampler = lambda x, y: y + resampler = lambda x, y: y # noqa else: from pytential.symbolic.mappers import stringify_where raise ValueError('unknown `where` identifier: {}'.format( stringify_where(where))) - density = self.rec(expr.density) - return resampler(self.queue, density) + return resampler(self.queue, operand) # }}} @@ -441,7 +449,6 @@ class GeometryCollection(object): if isinstance(where, sym.QBXSourceQuadStage2): return lpot.quad_stage2_density_discr - from pytential.symbolic.mappers import stringify_where raise ValueError('unknown `where` identifier: {}'.format( sym.stringify_where(where))) @@ -551,7 +558,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 InterpolationPreprocessor expr = _prepare_expr(places, expr) + expr = InterpolationPreprocessor()(expr) return BoundExpression(places, expr) diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 79b641ba..a2479a16 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -81,7 +81,6 @@ class IdentityMapper(IdentityMapperBase): map_node_coordinate_component = map_ones map_parametrization_gradient = map_ones map_parametrization_derivative = map_ones - map_interpolation = map_ones # }}} @@ -104,19 +103,20 @@ class IdentityMapper(IdentityMapperBase): for name, arg_expr in expr.kernel_arguments.items() )) + def map_interpolation(self, expr): + return type(expr)(self.rec(expr.operand), expr.source, expr.target) + class CombineMapper(CombineMapperBase): def map_node_sum(self, expr): return self.rec(expr.operand) - def map_interpolation(self, expr): - return self.rec(expr.density) - map_node_max = map_node_sum map_num_reference_derivative = map_node_sum 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( @@ -293,7 +293,7 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper): if target is None: target = self.default_where - return type(expr)(expr.density, source, target) + return type(expr)(expr.operand, source, target) def operand_rec(self, expr): return self.rec(expr) @@ -409,12 +409,19 @@ class UnregularizedPreprocessor(IdentityMapper): # {{{ QBX preprocessor class InterpolationPreprocessor(IdentityMapper): - def mat_int_g(self, expr): + def map_int_g(self, expr): if isinstance(expr.source, prim.QBXSourceQuadStage2): return expr density = prim.Interpolation(self.rec(expr.density), expr.source) - return expr.copy(density=density) + kernel_arguments = dict( + (name, prim.Interpolation(self.rec(arg_expr), expr.source)) + for name, arg_expr in expr.kernel_arguments.items()) + + return expr.copy( + kernel=expr.kernel, + density=density, + kernel_arguments=kernel_arguments) class QBXPreprocessor(IdentityMapper): @@ -578,6 +585,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 a19d97c8..7100df0e 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -1018,19 +1018,28 @@ class _NoArgSentinel(object): class Interpolation(Expression): - """Interpolate density from *source* to *target* discretization.""" + """Interpolate quantity from *source* to *target* discretization.""" - init_arg_names = ("density", "source", "target") + init_arg_names = ("operand", "source", "target") - def __init__(self, density, source, target=None): - self.density = density + def __new__(cls, operand, source, target=None): + if isinstance(operand, np.ndarray): + def make_op(operand_i): + return cls(operand_i, source, target) + + return componentwise(make_op, operand) + else: + return Expression.__new__(cls) + + def __init__(self, operand, source, target=None): + self.operand = operand self.source = source self.target = target if self.target is None: self.target = QBXSourceQuadStage2(DEFAULT_SOURCE) def __getinitargs__(self): - return (self.density, self.source, self.target) + return (self.operand, self.source, self.target) mapper_method = intern("map_interpolation") diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py index 1e048b7f..b65b5263 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 33efe9f9..4b942f90 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -211,7 +211,6 @@ def test_interpolation(ctx_factory, where): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - ndim = 2 nelements = 32 target_order = 5 qbx_order = 4 @@ -241,7 +240,7 @@ def test_interpolation(ctx_factory, where): 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 = 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) -- GitLab From ded6487524cee15fa47ccac6b257f1cd8f747d9e Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Wed, 28 Nov 2018 21:21:56 -0600 Subject: [PATCH 04/11] execution: whitelist in map_interpolation --- pytential/symbolic/execution.py | 65 +++++++++++++++++--------------- pytential/symbolic/mappers.py | 24 ++++++++---- pytential/symbolic/primitives.py | 1 + 3 files changed, 51 insertions(+), 39 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 0d60ac29..650397bb 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -174,34 +174,37 @@ class EvaluationMapper(EvaluationMapperBase): return source.map_quad_kernel_op(expr, self.bound_expr, self.rec) def map_interpolation(self, expr): - assert isinstance(expr.target, sym.QBXSourceQuadStage2) + 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 not isinstance(source, LayerPotentialSourceBase): - return operand + if isinstance(source, LayerPotentialSourceBase): + where = expr.source + if not isinstance(where, sym._QBXSource): + where = sym.QBXSourceStage1(where) + + if isinstance(where, sym.QBXSourceStage1): + resampler = source.resampler + elif isinstance(where, sym.QBXSourceStage2): + resampler = source.refined_interp_to_ovsmp_quad_connection + elif isinstance(where, sym.QBXSourceQuadStage2): + resampler = lambda x, y: y # noqa + else: + from pytential.symbolic.mappers import stringify_where + raise ValueError('unknown `where` identifier: {}'.format( + stringify_where(where))) + else: + raise TypeError("source must be a `LayerPotentialSourceBase`") - if not isinstance(operand, cl.array.Array): + if isinstance(operand, cl.array.Array): + return resampler(self.queue, operand) + elif isinstance(operand, (int, float, complex, np.number)): return operand - - where = expr.source - if not isinstance(where, sym._QBXSource): - where = sym.QBXSourceStage1(where) - - if isinstance(where, sym.QBXSourceStage1): - resampler = source.resampler - elif isinstance(where, sym.QBXSourceStage2): - resampler = source.refined_interp_to_ovsmp_quad_connection - elif isinstance(where, sym.QBXSourceQuadStage2): - resampler = lambda x, y: y # noqa else: - from pytential.symbolic.mappers import stringify_where - raise ValueError('unknown `where` identifier: {}'.format( - stringify_where(where))) - - return resampler(self.queue, operand) + raise TypeError("cannot interpolate `{}`".format(type(operand))) # }}} @@ -417,7 +420,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): @@ -434,11 +437,7 @@ 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) @@ -464,14 +463,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] @@ -560,7 +563,7 @@ def bind(places, expr, auto_where=None): from pytential.symbolic.mappers import InterpolationPreprocessor expr = _prepare_expr(places, expr) - expr = InterpolationPreprocessor()(expr) + expr = InterpolationPreprocessor(places)(expr) return BoundExpression(places, expr) diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index a2479a16..b3473cc9 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -409,20 +409,28 @@ class UnregularizedPreprocessor(IdentityMapper): # {{{ QBX preprocessor class InterpolationPreprocessor(IdentityMapper): + def __init__(self, places): + self.places = places + def map_int_g(self, expr): if isinstance(expr.source, prim.QBXSourceQuadStage2): return expr - density = prim.Interpolation(self.rec(expr.density), expr.source) - kernel_arguments = dict( - (name, prim.Interpolation(self.rec(arg_expr), expr.source)) - for name, arg_expr in expr.kernel_arguments.items()) + from pytential.source import LayerPotentialSourceBase + source = self.places[expr.source] - return expr.copy( - kernel=expr.kernel, - density=density, - kernel_arguments=kernel_arguments) + if isinstance(source, LayerPotentialSourceBase): + density = prim.Interpolation(self.rec(expr.density), expr.source) + kernel_arguments = dict( + (name, prim.Interpolation(self.rec(arg_expr), expr.source)) + 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): diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 7100df0e..1b8fa1ff 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -114,6 +114,7 @@ Discretization properties ^^^^^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: QWeight +.. autoclass:: Interpolation .. autofunction:: nodes .. autofunction:: parametrization_derivative .. autofunction:: parametrization_derivative_matrix -- GitLab From d895aed0da1b9222a27a8a0575cf565beaf12590 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Thu, 29 Nov 2018 09:12:48 -0600 Subject: [PATCH 05/11] mappers: fix flake8 --- pytential/symbolic/mappers.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index b3473cc9..b2dbbfa6 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -423,7 +423,7 @@ class InterpolationPreprocessor(IdentityMapper): density = prim.Interpolation(self.rec(expr.density), expr.source) kernel_arguments = dict( (name, prim.Interpolation(self.rec(arg_expr), expr.source)) - for name, arg_expr in expr.kernel_arguments.items()) + for name, arg_expr in expr.kernel_arguments.items()) expr = expr.copy( kernel=expr.kernel, @@ -432,6 +432,7 @@ class InterpolationPreprocessor(IdentityMapper): return expr + class QBXPreprocessor(IdentityMapper): def __init__(self, source_name, places): self.source_name = source_name -- GitLab From 18f132d332ae2f0ded08569e885226f1af3d81de Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Thu, 29 Nov 2018 09:16:04 -0600 Subject: [PATCH 06/11] primitives: add interpolation docs to a new section --- pytential/symbolic/primitives.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 1b8fa1ff..bb98b2db 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -114,7 +114,6 @@ Discretization properties ^^^^^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: QWeight -.. autoclass:: Interpolation .. autofunction:: nodes .. autofunction:: parametrization_derivative .. autofunction:: parametrization_derivative_matrix @@ -142,6 +141,11 @@ Elementary numerics .. autofunction:: mean .. autoclass:: IterativeInverse +Operators +^^^^^^^^^ + +.. autoclass:: Interpolation + Geometric Calculus (based on Geometric/Clifford Algebra) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -- GitLab From a740b81b467e5f23056c911d6ff2cdada931cd86 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sat, 1 Dec 2018 19:43:48 -0600 Subject: [PATCH 07/11] primitives: change argument order of Interpolation --- pytential/symbolic/execution.py | 4 ++-- pytential/symbolic/mappers.py | 14 +++++++++----- pytential/symbolic/primitives.py | 14 ++++++-------- test/test_symbolic.py | 17 +++++++++-------- 4 files changed, 26 insertions(+), 23 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 650397bb..0de7aecc 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -561,9 +561,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 InterpolationPreprocessor + from pytential.symbolic.mappers import QBXInterpolationPreprocessor expr = _prepare_expr(places, expr) - expr = InterpolationPreprocessor(places)(expr) + expr = QBXInterpolationPreprocessor(places)(expr) return BoundExpression(places, expr) diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index b2dbbfa6..bfcc551e 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -104,7 +104,7 @@ class IdentityMapper(IdentityMapperBase): )) def map_interpolation(self, expr): - return type(expr)(self.rec(expr.operand), expr.source, expr.target) + return type(expr)(expr.source, expr.target, self.rec(expr.operand)) class CombineMapper(CombineMapperBase): @@ -293,7 +293,7 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper): if target is None: target = self.default_where - return type(expr)(expr.operand, source, target) + return type(expr)(source, target, expr.operand) def operand_rec(self, expr): return self.rec(expr) @@ -408,7 +408,7 @@ class UnregularizedPreprocessor(IdentityMapper): # {{{ QBX preprocessor -class InterpolationPreprocessor(IdentityMapper): +class QBXInterpolationPreprocessor(IdentityMapper): def __init__(self, places): self.places = places @@ -420,9 +420,13 @@ class InterpolationPreprocessor(IdentityMapper): source = self.places[expr.source] if isinstance(source, LayerPotentialSourceBase): - density = prim.Interpolation(self.rec(expr.density), expr.source) + target = prim.QBXSourceQuadStage2(prim.DEFAULT_SOURCE) + + density = prim.Interpolation( + expr.source, target, self.rec(expr.density)) kernel_arguments = dict( - (name, prim.Interpolation(self.rec(arg_expr), expr.source)) + (name, prim.Interpolation( + expr.source, target, self.rec(arg_expr))) for name, arg_expr in expr.kernel_arguments.items()) expr = expr.copy( diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index bb98b2db..ec129369 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -1025,26 +1025,24 @@ class _NoArgSentinel(object): class Interpolation(Expression): """Interpolate quantity from *source* to *target* discretization.""" - init_arg_names = ("operand", "source", "target") + init_arg_names = ("source", "target", "operand") - def __new__(cls, operand, source, target=None): + def __new__(cls, source, target, operand): if isinstance(operand, np.ndarray): def make_op(operand_i): - return cls(operand_i, source, target) + return cls(source, target, operand_i) return componentwise(make_op, operand) else: return Expression.__new__(cls) - def __init__(self, operand, source, target=None): - self.operand = operand + def __init__(self, source, target, operand): self.source = source self.target = target - if self.target is None: - self.target = QBXSourceQuadStage2(DEFAULT_SOURCE) + self.operand = operand def __getinitargs__(self): - return (self.operand, self.source, self.target) + return (self.source, self.target, self.operand) mapper_method = intern("map_interpolation") diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 4b942f90..d709f7f0 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -201,18 +201,18 @@ def test_expr_pickling(): # }}} -@pytest.mark.parametrize("where", [ +@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, where): +def test_interpolation(ctx_factory, source): ctx = ctx_factory() queue = cl.CommandQueue(ctx) nelements = 32 - target_order = 5 + target_order = 7 qbx_order = 4 mesh = make_curve_mesh(starfish, @@ -227,14 +227,15 @@ def test_interpolation(ctx_factory, where): 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(sigma_sym, where)) - bound_op = bind(qbx, op_sym, auto_where=(where, sym.DEFAULT_TARGET)) + 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(where, sym.QBXSourceStage2): + if isinstance(source, sym.QBXSourceStage2): nodes = qbx.stage2_density_discr.nodes().get(queue) - elif isinstance(where, sym.QBXSourceQuadStage2): + elif isinstance(source, sym.QBXSourceQuadStage2): nodes = quad2_nodes else: nodes = qbx.density_discr.nodes().get(queue) @@ -244,7 +245,7 @@ def test_interpolation(ctx_factory, where): sigma_quad2_interp = bound_op(queue, sigma=sigma_dev).get(queue) error = la.norm(sigma_quad2_interp - sigma_quad2) / la.norm(sigma_quad2) - print(error) + assert error < 1.0e-10 # You can test individual routines by typing -- GitLab From 8474b7cbe766c95afd20aa60d38963ab24d5f2ba Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 21 May 2019 07:51:19 -0500 Subject: [PATCH 08/11] execution: fix lack of import in case of error --- pytential/symbolic/execution.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 0de7aecc..901db436 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -448,8 +448,9 @@ class GeometryCollection(object): if isinstance(where, sym.QBXSourceQuadStage2): return lpot.quad_stage2_density_discr + from pytential.symbolic.mappers import stringify_where raise ValueError('unknown `where` identifier: {}'.format( - sym.stringify_where(where))) + stringify_where(where))) def get_discretization(self, where): """ -- GitLab From 3ee87fbc5065945a7317c62e79c2e54c2f706d15 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 21 May 2019 23:45:41 +0200 Subject: [PATCH 09/11] Apply suggestion to pytential/symbolic/execution.py --- pytential/symbolic/execution.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 901db436..f5a29c35 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -182,20 +182,22 @@ class EvaluationMapper(EvaluationMapperBase): from pytential.source import LayerPotentialSourceBase if isinstance(source, LayerPotentialSourceBase): - where = expr.source - if not isinstance(where, sym._QBXSource): - where = sym.QBXSourceStage1(where) + sym_source = expr.source + if not isinstance(sym_source, sym._QBXSource): + sym_source = sym.QBXSourceStage1(sym_source) - if isinstance(where, sym.QBXSourceStage1): + if isinstance(sym_source, sym.QBXSourceStage1): resampler = source.resampler - elif isinstance(where, sym.QBXSourceStage2): + elif isinstance(sym_source, sym.QBXSourceStage2): resampler = source.refined_interp_to_ovsmp_quad_connection - elif isinstance(where, sym.QBXSourceQuadStage2): + 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: {}'.format( - stringify_where(where))) + raise ValueError( + "unknown `where` identifier in " + "interpolation source: {}".format( + stringify_where(sym_source))) else: raise TypeError("source must be a `LayerPotentialSourceBase`") -- GitLab From 62bf2586707cb9f6e3e9942920e5de1d5f9dca70 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 22 May 2019 00:13:39 +0200 Subject: [PATCH 10/11] Apply suggestion to pytential/symbolic/mappers.py --- pytential/symbolic/mappers.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index bfcc551e..41545353 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -420,13 +420,13 @@ class QBXInterpolationPreprocessor(IdentityMapper): source = self.places[expr.source] if isinstance(source, LayerPotentialSourceBase): - target = prim.QBXSourceQuadStage2(prim.DEFAULT_SOURCE) + stage2_source = prim.QBXSourceQuadStage2(prim.DEFAULT_SOURCE) density = prim.Interpolation( - expr.source, target, self.rec(expr.density)) + expr.source, stage2_source, self.rec(expr.density)) kernel_arguments = dict( (name, prim.Interpolation( - expr.source, target, self.rec(arg_expr))) + expr.source, stage2_source, self.rec(arg_expr))) for name, arg_expr in expr.kernel_arguments.items()) expr = expr.copy( -- GitLab From c6bd5b43889d572b28ef70ebcb0de5d5649d6747 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Wed, 22 May 2019 01:39:03 +0200 Subject: [PATCH 11/11] Placate flake8 --- pytential/symbolic/execution.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index f5a29c35..9a9bcb05 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -197,7 +197,7 @@ class EvaluationMapper(EvaluationMapperBase): raise ValueError( "unknown `where` identifier in " "interpolation source: {}".format( - stringify_where(sym_source))) + stringify_where(sym_source))) else: raise TypeError("source must be a `LayerPotentialSourceBase`") -- GitLab