From 3d0d63675f86a997053dbc1224948e6d388b22d1 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 2 Jun 2021 21:06:53 -0500 Subject: [PATCH 1/5] use geometry collection to evaluate matrix intg arguments That's where all the geometry is cached, so we might as well reuse it. --- pytential/symbolic/matrix.py | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 83bb8ed2e..f77ba9d09 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -40,21 +40,25 @@ def is_zero(x): return isinstance(x, (int, float, complex, np.number)) and x == 0 -def _get_layer_potential_args(mapper, expr, include_args=None): +def _get_layer_potential_args(actx, places, expr, context=None, include_args=None): """ :arg mapper: a :class:`~pytential.symbolic.matrix.MatrixBuilderBase`. :arg expr: symbolic layer potential expression. :return: a mapping of kernel arguments evaluated by the *mapper*. """ + from pytential import bind + if context is None: + context = {} kernel_args = {} for arg_name, arg_expr in expr.kernel_arguments.items(): - if (include_args is not None - and arg_name not in include_args): + if include_args is not None and arg_name not in include_args: continue - kernel_args[arg_name] = mapper.rec(arg_expr) + kernel_args[arg_name] = flatten( + bind(places, arg_expr)(actx, **context), + strict=False) return kernel_args @@ -376,7 +380,8 @@ def map_int_g(self, expr): raise NotImplementedError("layer potentials on non-variables") actx = self.array_context - kernel_args = _get_layer_potential_args(self, expr) + kernel_args = _get_layer_potential_args( + actx, self.places, expr, context=self.context) local_expn = lpot_source.get_expansion_for_qbx_direct_eval( kernel.get_base_kernel(), (expr.target_kernel,)) @@ -448,8 +453,9 @@ def map_int_g(self, expr): kernel_args = {arg.loopy_arg.name for arg in kernel_args} actx = self.array_context - kernel_args = _get_layer_potential_args(self, - expr, include_args=kernel_args) + kernel_args = _get_layer_potential_args( + actx, self.places, expr, context=self.context, + include_args=kernel_args) if self.exclude_self: kernel_args["target_to_source"] = actx.from_numpy( np.arange(0, target_discr.ndofs, dtype=np.int64) @@ -510,7 +516,8 @@ def map_int_g(self, expr): raise NotImplementedError actx = self.array_context - kernel_args = _get_layer_potential_args(self._mat_mapper, expr) + kernel_args = _get_layer_potential_args( + actx, self.places, expr, context=self.context) local_expn = lpot_source.get_expansion_for_qbx_direct_eval( kernel.get_base_kernel(), (expr.target_kernel,)) @@ -585,8 +592,9 @@ def map_int_g(self, expr): kernel_args = {arg.loopy_arg.name for arg in kernel_args} actx = self.array_context - kernel_args = _get_layer_potential_args(self._mat_mapper, - expr, include_args=kernel_args) + kernel_args = _get_layer_potential_args( + actx, self.places, expr, context=self.context, + include_args=kernel_args) if self.exclude_self: kernel_args["target_to_source"] = actx.from_numpy( np.arange(0, target_discr.ndofs, dtype=np.int64) From ba9729228ae1f731887b3469f7c73f36c42cbef3 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 2 Jun 2021 21:10:03 -0500 Subject: [PATCH 2/5] mark matrix p2p flatten with strict=False --- pytential/symbolic/matrix.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index f77ba9d09..dacf5a67b 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -401,7 +401,7 @@ def map_int_g(self, expr): dofdesc=expr.target))(actx) _, (mat,) = mat_gen(actx.queue, - targets=flatten(thaw(target_discr.nodes(), actx)), + targets=flatten(thaw(target_discr.nodes(), actx), strict=False), sources=flatten(thaw(source_discr.nodes(), actx)), centers=flatten(centers), expansion_radii=flatten(radii), @@ -411,7 +411,7 @@ def map_int_g(self, expr): waa = bind(self.places, sym.weights_and_area_elements( source_discr.ambient_dim, dofdesc=expr.source))(actx) - mat[:, :] *= actx.to_numpy(flatten(waa)) + mat[:, :] *= flatten_to_numpy(actx, waa) mat = mat.dot(rec_density) result += mat @@ -468,8 +468,8 @@ def map_int_g(self, expr): exclude_self=self.exclude_self) _, (mat,) = mat_gen(actx.queue, - targets=flatten(thaw(target_discr.nodes(), actx)), - sources=flatten(thaw(source_discr.nodes(), actx)), + targets=flatten(thaw(target_discr.nodes(), actx), strict=False), + sources=flatten(thaw(source_discr.nodes(), actx), strict=False), **kernel_args) result += actx.to_numpy(mat).dot(rec_density) @@ -536,7 +536,7 @@ def map_int_g(self, expr): dofdesc=expr.target))(actx) _, (mat,) = mat_gen(actx.queue, - targets=flatten(thaw(target_discr.nodes(), actx)), + targets=flatten(thaw(target_discr.nodes(), actx), strict=False), sources=flatten(thaw(source_discr.nodes(), actx)), centers=flatten(centers), expansion_radii=flatten(radii), @@ -607,8 +607,8 @@ def map_int_g(self, expr): exclude_self=self.exclude_self) _, (mat,) = mat_gen(actx.queue, - targets=flatten(thaw(target_discr.nodes(), actx)), - sources=flatten(thaw(source_discr.nodes(), actx)), + targets=flatten(thaw(target_discr.nodes(), actx), strict=False), + sources=flatten(thaw(source_discr.nodes(), actx), strict=False), index_set=self.index_set, **kernel_args) From 4b37a2639bb61b8f702e95c84c71693f21e763c4 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 2 Jun 2021 21:11:58 -0500 Subject: [PATCH 3/5] use base_kernel when evaluating p2p --- pytential/symbolic/matrix.py | 28 ++++++++-------------------- 1 file changed, 8 insertions(+), 20 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index dacf5a67b..070e2fc74 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -259,20 +259,6 @@ def __init__(self, actx, dep_expr, other_dep_exprs, places, context) self.index_set = index_set - @property - @memoize_method - def _mat_mapper(self): - # mat_mapper is used to compute any kernel arguments that needs to - # be computed on the full discretization, ignoring our index_set, - # e.g the normal in a double layer potential - - return MatrixBuilderBase(self.array_context, - self.dep_expr, - self.other_dep_exprs, - self.dep_source, - self.dep_discr, - self.places, self.context) - @property @memoize_method def _blk_mapper(self): @@ -449,7 +435,8 @@ def map_int_g(self, expr): raise NotImplementedError("layer potentials on non-variables") # NOTE: copied from pytential.symbolic.primitives.IntG - kernel_args = kernel.get_args() + kernel.get_source_args() + base_kernel = kernel.get_base_kernel() + kernel_args = base_kernel.get_args() + base_kernel.get_source_args() kernel_args = {arg.loopy_arg.name for arg in kernel_args} actx = self.array_context @@ -463,8 +450,8 @@ def map_int_g(self, expr): from sumpy.p2p import P2PMatrixGenerator mat_gen = P2PMatrixGenerator(actx.context, - source_kernels=(kernel,), - target_kernels=(expr.target_kernel,), + source_kernels=(base_kernel,), + target_kernels=(expr.target_kernel.get_base_kernel(),), exclude_self=self.exclude_self) _, (mat,) = mat_gen(actx.queue, @@ -588,7 +575,8 @@ def map_int_g(self, expr): raise NotImplementedError # NOTE: copied from pytential.symbolic.primitives.IntG - kernel_args = kernel.get_args() + kernel.get_source_args() + base_kernel = kernel.get_base_kernel() + kernel_args = base_kernel.get_args() + base_kernel.get_source_args() kernel_args = {arg.loopy_arg.name for arg in kernel_args} actx = self.array_context @@ -602,8 +590,8 @@ def map_int_g(self, expr): from sumpy.p2p import P2PMatrixBlockGenerator mat_gen = P2PMatrixBlockGenerator(actx.context, - source_kernels=(kernel,), - target_kernels=(expr.target_kernel,), + source_kernels=(base_kernel,), + target_kernels=(expr.target_kernel.get_base_kernel(),), exclude_self=self.exclude_self) _, (mat,) = mat_gen(actx.queue, From 9bc58ec712fc7eb9af6754f71a35ca032f2b1d26 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 2 Jun 2021 21:19:19 -0500 Subject: [PATCH 4/5] allow weighting p2p matrices by area elements --- pytential/symbolic/matrix.py | 39 ++++++++++++++++++++++++++++-------- 1 file changed, 31 insertions(+), 8 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 070e2fc74..79dc93172 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -398,9 +398,8 @@ def map_int_g(self, expr): source_discr.ambient_dim, dofdesc=expr.source))(actx) mat[:, :] *= flatten_to_numpy(actx, waa) - mat = mat.dot(rec_density) - result += mat + result += mat @ rec_density return result @@ -411,12 +410,14 @@ def map_int_g(self, expr): class P2PMatrixBuilder(MatrixBuilderBase): def __init__(self, actx, dep_expr, other_dep_exprs, - dep_source, dep_discr, places, context, exclude_self=True): + dep_source, dep_discr, places, context, + weighted=False, exclude_self=True): super().__init__(actx, dep_expr, other_dep_exprs, dep_source, dep_discr, places, context) self.exclude_self = exclude_self + self.weighted = weighted def map_int_g(self, expr): source_discr = self.places.get_discretization( @@ -458,8 +459,18 @@ def map_int_g(self, expr): targets=flatten(thaw(target_discr.nodes(), actx), strict=False), sources=flatten(thaw(source_discr.nodes(), actx), strict=False), **kernel_args) + mat = actx.to_numpy(mat) + + from meshmode.discretization import Discretization + if self.weighted and isinstance(source_discr, Discretization): + from pytential import bind, sym + waa = bind(self.places, sym.weights_and_area_elements( + source_discr.ambient_dim, + dofdesc=expr.source))(actx) - result += actx.to_numpy(mat).dot(rec_density) + mat[:, :] *= flatten_to_numpy(actx, waa) + + result += mat @ rec_density return result @@ -493,7 +504,6 @@ def map_int_g(self, expr): raise NotImplementedError result = 0 - for kernel, density in zip(expr.source_kernels, expr.densities): rec_density = self._blk_mapper.rec(density) if is_zero(rec_density): @@ -536,18 +546,21 @@ def map_int_g(self, expr): waa = flatten(waa) mat *= waa[self.index_set.linear_col_indices] - result += rec_density * actx.to_numpy(mat) + result += actx.to_numpy(mat) * rec_density return result class FarFieldBlockBuilder(MatrixBlockBuilderBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, - places, index_set, context, exclude_self=False): + places, index_set, context, + weighted=False, exclude_self=False): super().__init__(queue, dep_expr, other_dep_exprs, dep_source, dep_discr, places, index_set, context) + self.exclude_self = exclude_self + self.weighted = weighted def get_dep_variable(self): queue = self.array_context.queue @@ -600,7 +613,17 @@ def map_int_g(self, expr): index_set=self.index_set, **kernel_args) - result += rec_density * actx.to_numpy(mat) + from meshmode.discretization import Discretization + if self.weighted and isinstance(source_discr, Discretization): + from pytential import bind, sym + waa = bind(self.places, sym.weights_and_area_elements( + source_discr.ambient_dim, + dofdesc=expr.source))(actx) + waa = flatten(waa) + + mat *= waa[self.index_set.linear_col_indices] + + result += actx.to_numpy(mat) * rec_density return result From 48e6d1e3eab01d8a7c57be591dcfd9391fde4006 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 2 Jun 2021 21:22:48 -0500 Subject: [PATCH 5/5] allow non-self evaluation in block p2p --- pytential/symbolic/matrix.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 79dc93172..855b093a4 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -575,9 +575,6 @@ def map_int_g(self, expr): target_discr = self.places.get_discretization( expr.target.geometry, expr.target.discr_stage) - if source_discr is not target_discr: - raise NotImplementedError - result = 0 for kernel, density in zip(expr.source_kernels, expr.densities): rec_density = self._blk_mapper.rec(density)