From e27a43cf4cd88fc8a223c25b190e73eeece434d0 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 24 Jul 2019 09:30:40 -0500 Subject: [PATCH 1/3] add mapper to upsample when differentiating on quad_stage2. The mapper looks for a num_reference_derivative on `quad_stage2` and replaces it by one on `stage2` wrapped inside an interpolation to `quad_stage2`. --- pytential/qbx/__init__.py | 9 +-- pytential/symbolic/execution.py | 11 ++-- pytential/symbolic/mappers.py | 59 +++++++++++++++--- pytential/symbolic/matrix.py | 102 +++++++++++++------------------ pytential/symbolic/primitives.py | 10 +-- test/test_matrix.py | 24 +++++--- 6 files changed, 119 insertions(+), 96 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index a29d20bd..c6479606 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -350,11 +350,12 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): def weights_and_area_elements(self): from pytential import bind, sym with cl.CommandQueue(self.cl_context) as queue: - waa = bind(self, sym.weights_and_area_elements( - self.ambient_dim, - where=sym.QBX_SOURCE_QUAD_STAGE2))(queue) + where = sym.DOFDescriptor( + sym.DEFAULT_SOURCE, + discr=sym.QBX_SOURCE_QUAD_STAGE2) - return waa.with_queue(None) + return bind(self, sym.weights_and_area_elements( + self.ambient_dim, where=where))(queue).with_queue(None) # }}} diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index 1a6c8873..0ed5d4d9 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -322,7 +322,7 @@ def _prepare_domains(nresults, places, domains, default_domain): return domains -def _prepare_expr(places, expr): +def _prepare_expr(places, expr, interpolate=True): """ :arg places: :class:`pytential.symbolic.execution.GeometryCollection`. :arg expr: a symbolic expression. @@ -332,7 +332,9 @@ def _prepare_expr(places, expr): from pytential.source import LayerPotentialSourceBase from pytential.symbolic.mappers import ( - ToTargetTagger, DerivativeBinder) + ToTargetTagger, + DerivativeBinder, + InterpolationPreprocessor) expr = ToTargetTagger(*places.auto_where)(expr) expr = DerivativeBinder()(expr) @@ -341,6 +343,8 @@ def _prepare_expr(places, expr): if isinstance(place, LayerPotentialSourceBase): expr = place.preprocess_optemplate(name, places, expr) + if interpolate: + expr = InterpolationPreprocessor(places)(expr) return expr # }}} @@ -566,10 +570,7 @@ 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 426073b7..8ada3846 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -219,12 +219,7 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper): else: return expr - def map_q_weight(self, expr): - dd = prim.as_dofdesc(expr.where) - if dd.where is None: - return type(expr)(where=dd.copy(where=self.default_source)) - else: - return expr + map_q_weight = map_ones def map_parametrization_derivative_component(self, expr): dd = prim.as_dofdesc(expr.where) @@ -423,17 +418,59 @@ class UnregularizedPreprocessor(IdentityMapper): # {{{ QBX preprocessor -class QBXInterpolationPreprocessor(IdentityMapper): - def __init__(self, places): +class InterpolationDiscretizationTagger(IdentityMapper): + def __init__(self, discr): + self.discr = discr + + def map_node_coordinate_component(self, expr): + dd = prim.as_dofdesc(expr.where) + if dd.discr == self.discr: + return expr + else: + return type(expr)( + expr.ambient_axis, + dd.copy(discr=self.discr)) + + def map_num_reference_derivative(self, expr): + dd = prim.as_dofdesc(expr.where) + if dd.discr == self.discr: + return expr + else: + return type(expr)( + expr.ref_axes, + self.rec(expr.operand), + dd.copy(discr=self.discr)) + + +class InterpolationPreprocessor(IdentityMapper): + def __init__(self, places, discr=None): self.places = places + self.discr = discr or prim.QBX_SOURCE_STAGE2 + self.tagger = InterpolationDiscretizationTagger(self.discr) + + def map_num_reference_derivative(self, expr): + target_dd = prim.as_dofdesc(expr.where) + if target_dd.discr != prim.QBX_SOURCE_QUAD_STAGE2: + return expr + + from pytential.qbx import QBXLayerPotentialSource + lpot_source = self.places[target_dd] + if not isinstance(lpot_source, QBXLayerPotentialSource): + return expr + + source_dd = target_dd.copy(discr=self.discr) + return prim.Interpolation( + source_dd, + target_dd, + self.rec(self.tagger(expr))) def map_int_g(self, expr): source_dd = prim.as_dofdesc(expr.source) if source_dd.discr == prim.QBX_SOURCE_QUAD_STAGE2: return expr - lpot_source = self.places[expr.source] from pytential.qbx import QBXLayerPotentialSource + lpot_source = self.places[source_dd] if not isinstance(lpot_source, QBXLayerPotentialSource): return expr @@ -448,7 +485,9 @@ class QBXInterpolationPreprocessor(IdentityMapper): return expr.copy( kernel=expr.kernel, density=density, - kernel_arguments=kernel_arguments) + kernel_arguments=kernel_arguments, + source=target_dd, + target=expr.target) class QBXPreprocessor(IdentityMapper): diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 5b60fb31..f69ed44a 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -43,34 +43,6 @@ def is_zero(x): return isinstance(x, (int, float, complex, np.number)) and x == 0 -def _resample_arg(queue, source, x): - """ - :arg queue: a :class:`pyopencl.CommandQueue`. - :arg source: a :class:`pytential.source.LayerPotentialSourceBase` subclass. - If it is not a layer potential source, no resampling is done. - :arg x: a :class:`numpy.ndarray`. - - :return: a resampled :class:`numpy.ndarray` (see - :method:`pytential.source.LayerPotentialSourceBase.resampler`). - """ - - from pytential.source import LayerPotentialSourceBase - if not isinstance(source, LayerPotentialSourceBase): - return x - - if not isinstance(x, np.ndarray): - return x - - if len(x.shape) >= 2: - raise RuntimeError("matrix variables in kernel arguments") - - def resample(y): - return source.resampler(queue, cl.array.to_device(queue, y)).get(queue) - - from pytools.obj_array import with_object_array_or_scalar - return with_object_array_or_scalar(resample, x) - - def _get_layer_potential_args(mapper, expr, source): """ :arg mapper: a :class:`pytential.symbolic.matrix.MatrixBuilderBase`. @@ -82,8 +54,7 @@ def _get_layer_potential_args(mapper, expr, source): kernel_args = {} for arg_name, arg_expr in six.iteritems(expr.kernel_arguments): - rec_arg = mapper.rec(arg_expr) - kernel_args[arg_name] = _resample_arg(mapper.queue, source, rec_arg) + kernel_args[arg_name] = mapper.rec(arg_expr) return kernel_args @@ -106,9 +77,7 @@ def _get_kernel_args(mapper, kernel, expr, source): for arg_name, arg_expr in six.iteritems(expr.kernel_arguments): if arg_name not in inner_kernel_args: continue - - rec_arg = mapper.rec(arg_expr) - kernel_args[arg_name] = _resample_arg(mapper.queue, source, rec_arg) + kernel_args[arg_name] = mapper.rec(arg_expr) return kernel_args @@ -398,11 +367,32 @@ class MatrixBuilder(MatrixBuilderBase): super(MatrixBuilder, self).__init__(queue, dep_expr, other_dep_exprs, dep_source, dep_discr, places, context) + def map_interpolation(self, expr): + source_dd = sym.as_dofdesc(expr.source) + target_dd = sym.as_dofdesc(expr.target) + if target_dd.discr != sym.QBX_SOURCE_QUAD_STAGE2: + raise RuntimeError("can only interpolate to QBX_SOURCE_QUAD_STAGE2") + + operand = self.rec(expr.operand) + if isinstance(operand, (int, float, complex, np.number)): + return operand + elif isinstance(operand, np.ndarray) and operand.ndim == 1: + from pytential.symbolic.dof_connection import connection_from_dds + conn = connection_from_dds(self.places, + expr.source, expr.target) + + operand = cl.array.to_device(self.queue, operand) + return conn(self.queue, operand).get(self.queue) + elif isinstance(operand, np.ndarray) and operand.ndim == 2: + resampler = self.places[source_dd].direct_resampler + mat = resampler.full_resample_matrix(self.queue).get(self.queue) + return mat.dot(operand) + else: + raise RuntimeError('unknown operand type: {}'.format(type(operand))) + def map_int_g(self, expr): source_dd = sym.as_dofdesc(expr.source) target_dd = sym.as_dofdesc(expr.target) - if source_dd.discr is None: - source_dd = source_dd.copy(discr=sym.QBX_SOURCE_QUAD_STAGE2) lpot_source = self.places[source_dd] source_discr = self.places.get_discretization(source_dd) @@ -417,11 +407,7 @@ class MatrixBuilder(MatrixBuilderBase): raise NotImplementedError("layer potentials on non-variables") kernel = expr.kernel - if source_dd.discr == target_dd.discr: - # NOTE: passing None to avoid any resampling - kernel_args = _get_layer_potential_args(self, expr, None) - else: - kernel_args = _get_layer_potential_args(self, expr, lpot_source) + kernel_args = _get_layer_potential_args(self, expr, lpot_source) from sumpy.expansion.local import LineTaylorLocalExpansion local_expn = LineTaylorLocalExpansion(kernel, lpot_source.qbx_order) @@ -444,15 +430,6 @@ class MatrixBuilder(MatrixBuilderBase): waa = _get_weights_and_area_elements(self.queue, lpot_source, source_discr) mat[:, :] *= waa.get(self.queue) - - if source_dd.discr != target_dd.discr: - # NOTE: we only resample sources - assert target_discr.nnodes < source_discr.nnodes - - resampler = lpot_source.direct_resampler - resample_mat = resampler.full_resample_matrix(self.queue).get(self.queue) - mat = mat.dot(resample_mat) - mat = mat.dot(rec_density) return mat @@ -472,9 +449,12 @@ class P2PMatrixBuilder(MatrixBuilderBase): self.exclude_self = exclude_self def map_int_g(self, expr): - source = self.places[expr.source] - source_discr = self.places.get_discretization(expr.source) - target_discr = self.places.get_discretization(expr.target) + target_dd = sym.as_dofdesc(expr.target) + source_dd = sym.as_dofdesc(expr.source) + + source = self.places[source_dd] + source_discr = self.places.get_discretization(source_dd) + target_discr = self.places.get_discretization(target_dd) rec_density = self.rec(expr.density) if is_zero(rec_density): @@ -535,9 +515,12 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): - source = self.places[expr.source] - source_discr = self.places.get_discretization(expr.source) - target_discr = self.places.get_discretization(expr.target) + target_dd = sym.as_dofdesc(expr.target) + source_dd = sym.as_dofdesc(expr.source) + + source = self.places[source_dd] + source_discr = self.places.get_discretization(source_dd) + target_discr = self.places.get_discretization(target_dd) if source_discr is not target_discr: raise NotImplementedError() @@ -601,9 +584,12 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): return np.equal(tgtindices, srcindices).astype(np.float64) def map_int_g(self, expr): - source = self.places[expr.source] - source_discr = self.places.get_discretization(expr.source) - target_discr = self.places.get_discretization(expr.target) + target_dd = sym.as_dofdesc(expr.target) + source_dd = sym.as_dofdesc(expr.source) + + source = self.places[source_dd] + source_discr = self.places.get_discretization(source_dd) + target_discr = self.places.get_discretization(target_dd) if source_discr is not target_discr: raise NotImplementedError() diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 091c2925..d82a8d2f 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -964,16 +964,8 @@ def h_max(ambient_dim, dim=None, where=None): def weights_and_area_elements(ambient_dim, dim=None, where=None): where = as_dofdesc(where) - if where.discr == QBX_SOURCE_QUAD_STAGE2: - # quad_stage2_density_discr is not guaranteed to be usable for - # interpolation/differentiation. Use stage2_density_discr to find - # area elements instead, then upsample that. - source = where.copy(discr=QBX_SOURCE_STAGE2) - area = Interpolation(source, where, - area_element(ambient_dim, dim=dim, where=source)) - else: - area = area_element(ambient_dim, dim=dim, where=where) + area = area_element(ambient_dim, dim=dim, where=where) return cse(area * QWeight(where=where), "weights_area_elements", cse_scope.DISCRETIZATION) diff --git a/test/test_matrix.py b/test/test_matrix.py index 1800cec2..c667535c 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -286,7 +286,7 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, from pytential.symbolic.execution import GeometryCollection from pytential.symbolic.execution import _prepare_expr, _prepare_domains places = GeometryCollection(qbx) - expr = _prepare_expr(places, op) + expr = _prepare_expr(places, op, interpolate=False) domains = _prepare_domains(1, places, None, sym.DEFAULT_SOURCE) from pytential.symbolic.matrix import P2PMatrixBuilder @@ -362,7 +362,7 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id, from pytential.symbolic.execution import GeometryCollection, _prepare_expr places = GeometryCollection(qbx, auto_where=place_ids) - expr = _prepare_expr(places, op) + expr = _prepare_expr(places, op, interpolate=False) density_discr = places.get_discretization(place_ids[0]) index_set = _build_block_index(density_discr, factor=factor) @@ -441,17 +441,21 @@ def test_build_matrix_places(ctx_factory, place_ids, visualize=False): index_set = _build_block_index(source_discr, factor=0.6) - # build full QBX matrix - from pytential.symbolic.execution import build_matrix - qbx_mat = build_matrix(queue, qbx, op, u_sym, auto_where=place_ids) - qbx_mat = qbx_mat.get(queue) + from pytential.symbolic.execution import _prepare_expr + op = _prepare_expr(places, op, interpolate=False) - assert qbx_mat.shape == (target_discr.nnodes, source_discr.nnodes) + # build full QBX matrix + from pytential.symbolic.matrix import MatrixBuilder + mbuilder = MatrixBuilder(queue, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places[place_ids[0]], + dep_discr=places.get_discretization(place_ids[0]), + places=places, + context={}) + qbx_mat = mbuilder(op) # build full p2p matrix - from pytential.symbolic.execution import _prepare_expr - op = _prepare_expr(places, op) - from pytential.symbolic.matrix import P2PMatrixBuilder mbuilder = P2PMatrixBuilder(queue, dep_expr=u_sym, -- GitLab From 851411e17278930f5df5e7ec995d3ee68ae200d5 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 26 Jul 2019 20:06:06 -0500 Subject: [PATCH 2/3] small cleanups --- pytential/qbx/__init__.py | 7 ++----- pytential/symbolic/mappers.py | 14 +++++++------- pytential/symbolic/matrix.py | 2 +- 3 files changed, 10 insertions(+), 13 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 159b7b4b..f6975db2 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -357,12 +357,9 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): def weights_and_area_elements(self): from pytential import bind, sym with cl.CommandQueue(self.cl_context) as queue: - where = sym.DOFDescriptor( - sym.DEFAULT_SOURCE, - discr=sym.QBX_SOURCE_QUAD_STAGE2) - return bind(self, sym.weights_and_area_elements( - self.ambient_dim, where=where))(queue).with_queue(None) + self.ambient_dim, + where=sym.QBX_SOURCE_QUAD_STAGE2))(queue).with_queue(None) # }}} diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 21b64434..f1218ddb 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -416,7 +416,7 @@ class InterpolationDiscretizationTagger(IdentityMapper): def map_node_coordinate_component(self, expr): dd = prim.as_dofdesc(expr.where) - if dd.discr == self.discr: + if dd.discr is self.discr: return expr else: return type(expr)( @@ -425,7 +425,7 @@ class InterpolationDiscretizationTagger(IdentityMapper): def map_num_reference_derivative(self, expr): dd = prim.as_dofdesc(expr.where) - if dd.discr == self.discr: + if dd.discr is self.discr: return expr else: return type(expr)( @@ -435,14 +435,14 @@ class InterpolationDiscretizationTagger(IdentityMapper): class InterpolationPreprocessor(IdentityMapper): - def __init__(self, places, discr=None): + def __init__(self, places, from_discr=None): self.places = places - self.discr = discr or prim.QBX_SOURCE_STAGE2 - self.tagger = InterpolationDiscretizationTagger(self.discr) + self.from_discr = from_discr or prim.QBX_SOURCE_STAGE2 + self.tagger = InterpolationDiscretizationTagger(self.from_discr) def map_num_reference_derivative(self, expr): target_dd = prim.as_dofdesc(expr.where) - if target_dd.discr != prim.QBX_SOURCE_QUAD_STAGE2: + if target_dd.discr is not prim.QBX_SOURCE_QUAD_STAGE2: return expr from pytential.qbx import QBXLayerPotentialSource @@ -450,7 +450,7 @@ class InterpolationPreprocessor(IdentityMapper): if not isinstance(lpot_source, QBXLayerPotentialSource): return expr - source_dd = target_dd.copy(discr=self.discr) + source_dd = target_dd.copy(discr=self.from_discr) return prim.Interpolation( source_dd, target_dd, diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index f69ed44a..42f622f8 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -370,7 +370,7 @@ class MatrixBuilder(MatrixBuilderBase): def map_interpolation(self, expr): source_dd = sym.as_dofdesc(expr.source) target_dd = sym.as_dofdesc(expr.target) - if target_dd.discr != sym.QBX_SOURCE_QUAD_STAGE2: + if target_dd.discr is not sym.QBX_SOURCE_QUAD_STAGE2: raise RuntimeError("can only interpolate to QBX_SOURCE_QUAD_STAGE2") operand = self.rec(expr.operand) -- GitLab From caf31175eefe0e43bab18a5cd2354301cfb479e5 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 28 Jul 2019 13:28:31 -0500 Subject: [PATCH 3/3] docs and cleanup --- pytential/symbolic/execution.py | 7 ++- pytential/symbolic/mappers.py | 102 ++++++++++++++++++++----------- pytential/symbolic/matrix.py | 2 +- pytential/symbolic/primitives.py | 9 ++- test/test_matrix.py | 91 ++++++++++++++++----------- 5 files changed, 133 insertions(+), 78 deletions(-) diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index f263ce0f..f1f9fba1 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -372,7 +372,7 @@ def _prepare_domains(nresults, places, domains, default_domain): return domains -def _prepare_expr(places, expr, interpolate=True): +def _prepare_expr(places, expr): """ :arg places: :class:`pytential.symbolic.execution.GeometryCollection`. :arg expr: a symbolic expression. @@ -393,8 +393,9 @@ def _prepare_expr(places, expr, interpolate=True): if isinstance(place, LayerPotentialSourceBase): expr = place.preprocess_optemplate(name, places, expr) - if interpolate: - expr = InterpolationPreprocessor(places)(expr) + # NOTE: only insert interpolation operators after the layer potential + # operators were preprocessed to avoid any confusion + expr = InterpolationPreprocessor(places)(expr) return expr # }}} diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index d02f1a07..f7f00480 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -152,7 +152,7 @@ class DependencyMapper(DependencyMapperBase, Collector): class EvaluationMapper(EvaluationMapperBase): """Unlike :mod:`pymbolic.mapper.evaluation.EvaluationMapper`, this class does evaluation mostly to get :class:`pymbolic.geometric_algebra.MultiVector` - instances to to do their thing, and perhaps to automatically kill terms + instances to do their thing, and perhaps to automatically kill terms that are multiplied by zero. Otherwise it intends to largely preserve the structure of the input expression. """ @@ -199,7 +199,7 @@ class EvaluationMapper(EvaluationMapperBase): expr.scope) -# {{{ target/source tagging +# {{{ dofdesc tagging class LocationTagger(CSECachingMapperMixin, IdentityMapper): """Used internally by :class:`ToTargetTagger`. Tags all src/target taggable @@ -315,6 +315,46 @@ class ToTargetTagger(LocationTagger): self.operand_rec = LocationTagger(default_source, default_source=default_source) + +class DiscretizationStageTagger(IdentityMapper): + """Descends into an expression tree and changes the + :attr:`~pytential.symbolic.primitives.DOFDescriptor.discr_stage` to + :attr:`discr_stage`. + + .. attribute:: discr_stage + + The new discretization for the DOFs in the expression. For valid + values, see + :attr:`~pytential.symbolic.primitives.DOFDescriptor.discr_stage`. + """ + + def __init__(self, discr_stage): + if not (discr_stage == prim.QBX_SOURCE_STAGE1 + or discr_stage == prim.QBX_SOURCE_STAGE2 + or discr_stage == prim.QBX_SOURCE_QUAD_STAGE2): + raise ValueError('unknown discr stage tag: "{}"'.format(discr_stage)) + + self.discr_stage = discr_stage + + def map_node_coordinate_component(self, expr): + dofdesc = prim.as_dofdesc(expr.where) + if dofdesc.discr_stage == self.discr_stage: + return expr + + return type(expr)( + expr.ambient_axis, + dofdesc.copy(discr_stage=self.discr_stage)) + + def map_num_reference_derivative(self, expr): + dofdesc = prim.as_dofdesc(expr.where) + if dofdesc.discr_stage == self.discr_stage: + return expr + + return type(expr)( + expr.ref_axes, + self.rec(expr.operand), + dofdesc.copy(discr_stage=self.discr_stage)) + # }}} @@ -409,41 +449,29 @@ class UnregularizedPreprocessor(IdentityMapper): # }}} -# {{{ QBX preprocessor - -class InterpolationDiscretizationTagger(IdentityMapper): - def __init__(self, discr): - self.discr = discr - - def map_node_coordinate_component(self, expr): - dd = prim.as_dofdesc(expr.where) - if dd.discr is self.discr: - return expr - else: - return type(expr)( - expr.ambient_axis, - dd.copy(discr=self.discr)) - - def map_num_reference_derivative(self, expr): - dd = prim.as_dofdesc(expr.where) - if dd.discr is self.discr: - return expr - else: - return type(expr)( - expr.ref_axes, - self.rec(expr.operand), - dd.copy(discr=self.discr)) - +# {{{ interpolation preprocessor class InterpolationPreprocessor(IdentityMapper): - def __init__(self, places, from_discr=None): + """Handle expressions that require upsampling or downsampling by inserting + a :class:`~pytential.symbolic.primitives.Interpolation`. This is used to + + * do differentiation on + :attr:`~pytential.source.LayerPotentialSource.quad_stage2_density_discr`, + by performing it on + :attr:`~pytential.source.LayerPotentialSource.stage2_density_discr` and + upsampling. + * upsample layer potential sources to + :attr:`~pytential.source.LayerPotentialSource.quad_stage2_density_discr`, + """ + + def __init__(self, places): self.places = places - self.from_discr = from_discr or prim.QBX_SOURCE_STAGE2 - self.tagger = InterpolationDiscretizationTagger(self.from_discr) + self.from_discr_stage = prim.QBX_SOURCE_STAGE2 + self.tagger = DiscretizationStageTagger(self.from_discr_stage) def map_num_reference_derivative(self, expr): target_dd = prim.as_dofdesc(expr.where) - if target_dd.discr is not prim.QBX_SOURCE_QUAD_STAGE2: + if target_dd.discr_stage != prim.QBX_SOURCE_QUAD_STAGE2: return expr from pytential.qbx import QBXLayerPotentialSource @@ -451,7 +479,7 @@ class InterpolationPreprocessor(IdentityMapper): if not isinstance(lpot_source, QBXLayerPotentialSource): return expr - source_dd = target_dd.copy(discr=self.from_discr) + source_dd = target_dd.copy(discr_stage=self.from_discr_stage) return prim.Interpolation( source_dd, target_dd, @@ -459,11 +487,11 @@ class InterpolationPreprocessor(IdentityMapper): def map_int_g(self, expr): source = prim.as_dofdesc(expr.source) - if source.discr_stage == prim.QBX_SOURCE_QUAD_STAGE2: + if source.discr_stage is not None: return expr from pytential.qbx import QBXLayerPotentialSource - lpot_source = self.places[source_dd] + lpot_source = self.places[source] if not isinstance(lpot_source, QBXLayerPotentialSource): return expr @@ -479,9 +507,13 @@ class InterpolationPreprocessor(IdentityMapper): kernel=expr.kernel, density=density, kernel_arguments=kernel_arguments, - source=target_dd, + source=target, target=expr.target) +# }}} + + +# {{{ QBX preprocessor class QBXPreprocessor(IdentityMapper): def __init__(self, source_name, places): diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 42f622f8..a73e0c4d 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -370,7 +370,7 @@ class MatrixBuilder(MatrixBuilderBase): def map_interpolation(self, expr): source_dd = sym.as_dofdesc(expr.source) target_dd = sym.as_dofdesc(expr.target) - if target_dd.discr is not sym.QBX_SOURCE_QUAD_STAGE2: + if target_dd.discr_stage != sym.QBX_SOURCE_QUAD_STAGE2: raise RuntimeError("can only interpolate to QBX_SOURCE_QUAD_STAGE2") operand = self.rec(expr.operand) diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 6a0e4cbb..aa730beb 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -326,12 +326,11 @@ class DOFDescriptor(object): return not self.__eq__(other) def __repr__(self): + discr_stage = self.discr_stage \ + if self.discr_stage is None else self.discr_stage.__name__, + granularity = self.granularity.__name__ return '{}(geometry={}, stage={}, granularity={})'.format( - type(self).__name__, - self.geometry, - self.discr_stage \ - if self.discr_stage is None else self.discr_stage.__name__, - self.granularity.__name__) + type(self).__name__, self.geometry, discr_stage, granularity) def __str__(self): name = [] diff --git a/test/test_matrix.py b/test/test_matrix.py index 16b691c2..64b001bb 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -122,6 +122,8 @@ def _build_block_index(discr, nblks=10, factor=1.0): def _build_op(lpot_id, k=0, ndim=2, + source=sym.DEFAULT_SOURCE, + target=sym.DEFAULT_TARGET, qbx_forced_limit="avg"): from sumpy.kernel import LaplaceKernel, HelmholtzKernel @@ -132,7 +134,10 @@ def _build_op(lpot_id, knl = LaplaceKernel(ndim) knl_kwargs = {} - lpot_kwargs = {"qbx_forced_limit": qbx_forced_limit} + lpot_kwargs = { + "qbx_forced_limit": qbx_forced_limit, + "source": source, + "target": target} lpot_kwargs.update(knl_kwargs) if lpot_id == 1: # scalar single-layer potential @@ -278,23 +283,33 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, from sympy.core.cache import clear_cache clear_cache() + place_ids = ( + sym.DOFDescriptor( + geometry=sym.DEFAULT_SOURCE, + discr_stage=sym.QBX_SOURCE_STAGE1), + sym.DOFDescriptor( + geometry=sym.DEFAULT_TARGET, + discr_stage=sym.QBX_SOURCE_STAGE1), + ) target_order = 2 if ndim == 3 else 7 + qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) - op, u_sym, _ = _build_op(lpot_id, ndim=ndim) + op, u_sym, _ = _build_op(lpot_id, ndim=ndim, + source=place_ids[0], + target=place_ids[1]) index_set = _build_block_index(qbx.density_discr, factor=factor) from pytential.symbolic.execution import GeometryCollection - from pytential.symbolic.execution import _prepare_expr, _prepare_domains - places = GeometryCollection(qbx) - expr = _prepare_expr(places, op, interpolate=False) - domains = _prepare_domains(1, places, None, sym.DEFAULT_SOURCE) + from pytential.symbolic.execution import _prepare_expr + places = GeometryCollection(qbx, auto_where=place_ids) + expr = _prepare_expr(places, op) from pytential.symbolic.matrix import P2PMatrixBuilder mbuilder = P2PMatrixBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[domains[0]], - dep_discr=places.get_discretization(domains[0]), + dep_source=places[place_ids[0]], + dep_discr=places.get_discretization(place_ids[0]), places=places, context={}, exclude_self=True) @@ -304,8 +319,8 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id, mbuilder = FarFieldBlockBuilder(queue, dep_expr=u_sym, other_dep_exprs=[], - dep_source=places[domains[0]], - dep_discr=places.get_discretization(domains[0]), + dep_source=places[place_ids[0]], + dep_discr=places.get_discretization(place_ids[0]), places=places, index_set=index_set, context={}, @@ -347,24 +362,25 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id, from sympy.core.cache import clear_cache clear_cache() + place_ids = ( + sym.DOFDescriptor( + geometry=sym.DEFAULT_SOURCE, + discr_stage=sym.QBX_SOURCE_STAGE2), + sym.DOFDescriptor( + geometry=sym.DEFAULT_TARGET, + discr_stage=sym.QBX_SOURCE_STAGE2), + ) target_order = 2 if ndim == 3 else 7 + qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim) - op, u_sym, _ = _build_op(lpot_id, ndim=ndim) - - # NOTE: NearFieldBlockBuilder only does stage1/stage1 or stage2/stage2, - # so we need to hardcode the discr for MatrixBuilder too, since the - # defaults are different - source_dd = sym.DOFDescriptor( - geometry=sym.DEFAULT_SOURCE, - discr_stage=sym.QBX_SOURCE_STAGE2) - target_dd = sym.DOFDescriptor( - geometry=sym.DEFAULT_TARGET, - discr_stage=sym.QBX_SOURCE_STAGE2) - place_ids = (source_dd, target_dd) + op, u_sym, _ = _build_op(lpot_id, ndim=ndim, + source=place_ids[0], + target=place_ids[1], + qbx_forced_limit="avg") from pytential.symbolic.execution import GeometryCollection, _prepare_expr places = GeometryCollection(qbx, auto_where=place_ids) - expr = _prepare_expr(places, op, interpolate=False) + expr = _prepare_expr(places, op) density_discr = places.get_discretization(place_ids[0]) index_set = _build_block_index(density_discr, factor=factor) @@ -412,11 +428,11 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id, assert _max_block_error(mat, blk, index_set) < 1.0e-14 -@pytest.mark.parametrize('place_ids', - [(None, None), - (sym.QBX_SOURCE_STAGE1, sym.QBX_SOURCE_STAGE1), +@pytest.mark.parametrize(('source_discr_stage', 'target_discr_stage'), + [(sym.QBX_SOURCE_STAGE1, sym.QBX_SOURCE_STAGE1), (sym.QBX_SOURCE_QUAD_STAGE2, sym.QBX_SOURCE_QUAD_STAGE2)]) -def test_build_matrix_places(ctx_factory, place_ids, visualize=False): +def test_build_matrix_places(ctx_factory, + source_discr_stage, target_discr_stage, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) @@ -424,18 +440,25 @@ def test_build_matrix_places(ctx_factory, place_ids, visualize=False): from sympy.core.cache import clear_cache clear_cache() + qbx_forced_limit = -1 + place_ids = ( + sym.DOFDescriptor( + geometry=sym.DEFAULT_SOURCE, + discr_stage=source_discr_stage), + sym.DOFDescriptor( + geometry=sym.DEFAULT_TARGET, + discr_stage=target_discr_stage), + ) + + # build test operators qbx = _build_qbx_discr(queue, nelements=8, target_order=2, ndim=2, curve_f=partial(ellipse, 1.0)) - qbx_forced_limit = -1 op, u_sym, _ = _build_op(lpot_id=1, ndim=2, + source=place_ids[0], + target=place_ids[1], qbx_forced_limit=qbx_forced_limit) - place_ids = ( - sym.as_dofdesc(place_ids[0]).copy(geometry=sym.DEFAULT_SOURCE), - sym.as_dofdesc(place_ids[1]).copy(geometry=sym.DEFAULT_TARGET) - ) - from pytential.symbolic.execution import GeometryCollection places = GeometryCollection(qbx, auto_where=place_ids) source_discr = places.get_discretization(place_ids[0]) @@ -444,7 +467,7 @@ def test_build_matrix_places(ctx_factory, place_ids, visualize=False): index_set = _build_block_index(source_discr, factor=0.6) from pytential.symbolic.execution import _prepare_expr - op = _prepare_expr(places, op, interpolate=False) + op = _prepare_expr(places, op) # build full QBX matrix from pytential.symbolic.matrix import MatrixBuilder -- GitLab