diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 120de25239138d1663fa4ab2e51a7159bfd04624..f6975db2c4bd4f82494fd0c11199bb615fcc5f46 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -357,11 +357,9 @@ 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( + return bind(self, sym.weights_and_area_elements( self.ambient_dim, - where=sym.QBX_SOURCE_QUAD_STAGE2))(queue) - - return waa.with_queue(None) + where=sym.QBX_SOURCE_QUAD_STAGE2))(queue).with_queue(None) # }}} diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py index acbbc50e1d548d622727b023bdee40e70efaff23..f1f9fba13332df07925adb001c929e457a5da7cd 100644 --- a/pytential/symbolic/execution.py +++ b/pytential/symbolic/execution.py @@ -382,7 +382,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) @@ -391,6 +393,9 @@ def _prepare_expr(places, expr): if isinstance(place, LayerPotentialSourceBase): expr = place.preprocess_optemplate(name, places, expr) + # NOTE: only insert interpolation operators after the layer potential + # operators were preprocessed to avoid any confusion + expr = InterpolationPreprocessor(places)(expr) return expr # }}} @@ -618,10 +623,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 3bc719d9ede9c37f2c902eac1f159cb80ef2f582..f7f004806890cf2c098e4fabe9fd706b1e92bf32 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,19 +449,49 @@ class UnregularizedPreprocessor(IdentityMapper): # }}} -# {{{ QBX preprocessor +# {{{ interpolation preprocessor + +class InterpolationPreprocessor(IdentityMapper): + """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`, + """ -class QBXInterpolationPreprocessor(IdentityMapper): def __init__(self, places): self.places = places + 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_stage != 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_stage=self.from_discr_stage) + return prim.Interpolation( + source_dd, + target_dd, + self.rec(self.tagger(expr))) 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 - lpot_source = self.places[source] from pytential.qbx import QBXLayerPotentialSource + lpot_source = self.places[source] if not isinstance(lpot_source, QBXLayerPotentialSource): return expr @@ -436,9 +506,15 @@ class QBXInterpolationPreprocessor(IdentityMapper): return expr.copy( kernel=expr.kernel, density=density, - kernel_arguments=kernel_arguments) + kernel_arguments=kernel_arguments, + source=target, + target=expr.target) + +# }}} +# {{{ QBX preprocessor + class QBXPreprocessor(IdentityMapper): def __init__(self, source_name, places): self.source_name = source_name diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index a970d5451e28dfcad57a8f0c64803ea4b180d95c..a73e0c4d27f75c2bbf19bf2ae93439309f4cde7b 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_stage != 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_stage is None: - source_dd = source_dd.copy(discr_stage=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_stage == target_dd.discr_stage: - # 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_stage != target_dd.discr_stage: - # 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 72362ed21b9504d04707769392b71becd9015530..aa730beb47cc6b320366ce087177f57e580ade4a 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 = [] @@ -990,17 +989,7 @@ 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_stage == 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_stage=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 9fe20fa2aca49d86f39058c613a1f58dfeb0cfab..64b001bb35bf2fbbf825469321c3ec3543429560 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) + from pytential.symbolic.execution import _prepare_expr + places = GeometryCollection(qbx, auto_where=place_ids) expr = _prepare_expr(places, op) - domains = _prepare_domains(1, places, None, sym.DEFAULT_SOURCE) 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,20 +362,21 @@ 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) @@ -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]) @@ -443,17 +466,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) - - assert qbx_mat.shape == (target_discr.nnodes, source_discr.nnodes) - - # build full p2p matrix from pytential.symbolic.execution import _prepare_expr op = _prepare_expr(places, op) + # 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.matrix import P2PMatrixBuilder mbuilder = P2PMatrixBuilder(queue, dep_expr=u_sym,