From e7371f9f0d77fc48c75105e007fbbd2b07fa8ba8 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 22 Sep 2022 19:30:02 +0300 Subject: [PATCH] use the full double-layer in source-proxy skeletonization --- pytential/linalg/direct_solver_symbolic.py | 8 ++++++-- pytential/linalg/skeletonization.py | 15 ++++++++++----- pytential/symbolic/matrix.py | 18 ++++++------------ 3 files changed, 22 insertions(+), 19 deletions(-) diff --git a/pytential/linalg/direct_solver_symbolic.py b/pytential/linalg/direct_solver_symbolic.py index c57752840..5c11eee1d 100644 --- a/pytential/linalg/direct_solver_symbolic.py +++ b/pytential/linalg/direct_solver_symbolic.py @@ -54,12 +54,16 @@ def prepare_expr(places, exprs, auto_where=None): return _prepare_expr(places, exprs, auto_where=auto_where) -def prepare_proxy_expr(places, exprs, auto_where=None): +def prepare_proxy_expr( + places, exprs, + auto_where=None, + remove_transforms: bool = True): def _prepare_expr(expr): # remove all diagonal / non-operator terms in the expression expr = IntGTermCollector()(expr) # ensure all IntGs remove all the kernel derivatives - expr = KernelTransformationRemover()(expr) + if remove_transforms: + expr = KernelTransformationRemover()(expr) # ensure all IntGs have their source and targets set expr = DOFDescriptorReplacer(auto_where[0], auto_where[1])(expr) diff --git a/pytential/linalg/skeletonization.py b/pytential/linalg/skeletonization.py index 66d55c80a..88c2694f8 100644 --- a/pytential/linalg/skeletonization.py +++ b/pytential/linalg/skeletonization.py @@ -417,6 +417,7 @@ def make_skeletonization_wrangler( # internal _weighted_proxy: Optional[Union[bool, Tuple[bool, bool]]] = None, + _remove_source_transforms: bool = False, _proxy_source_cluster_builder: Optional[Callable[..., np.ndarray]] = None, _proxy_target_cluster_builder: Optional[Callable[..., np.ndarray]] = None, _neighbor_cluster_builder: Optional[Callable[..., np.ndarray]] = None, @@ -444,9 +445,13 @@ def make_skeletonization_wrangler( exprs = prepare_expr(places, exprs, auto_where) source_proxy_exprs = prepare_proxy_expr( - places, exprs, (auto_where[0], PROXY_SKELETONIZATION_TARGET)) + places, exprs, (auto_where[0], PROXY_SKELETONIZATION_TARGET), + remove_transforms=_remove_source_transforms) target_proxy_exprs = prepare_proxy_expr( - places, exprs, (PROXY_SKELETONIZATION_SOURCE, auto_where[1])) + places, exprs, (PROXY_SKELETONIZATION_SOURCE, auto_where[1]), + # NOTE: transforms are unconditionally removed here because the + # source would be the proxies, where we do not have normals, etc. + remove_transforms=True) # }}} @@ -456,7 +461,7 @@ def make_skeletonization_wrangler( weighted_sources = weighted_targets = True elif isinstance(_weighted_proxy, bool): weighted_sources = weighted_targets = _weighted_proxy - elif isinstance(_weighted_proxy, tuple): + elif isinstance(_weighted_proxy, tuple) and len(_weighted_proxy) == 2: weighted_sources, weighted_targets = _weighted_proxy else: raise ValueError(f"unknown value for weighting: '{_weighted_proxy}'") @@ -653,9 +658,9 @@ def _skeletonize_block_by_proxy_with_mats( if __debug__: isfinite = np.isfinite(tgt_mat) - assert np.all(isfinite), np.where(isfinite) + assert np.all(isfinite), np.where(~isfinite) isfinite = np.isfinite(src_mat) - assert np.all(isfinite), np.where(isfinite) + assert np.all(isfinite), np.where(~isfinite) # skeletonize target points k, idx, interp = interp_decomp(tgt_mat.T, rank=k, eps=id_eps) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index c24353dca..0484527bb 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -461,7 +461,6 @@ def map_int_g(self, expr): expr.target.geometry, expr.target.discr_stage) actx = self.array_context - target_base_kernel = expr.target_kernel.get_base_kernel() result = 0 for density, kernel in zip(expr.densities, expr.source_kernels): @@ -475,12 +474,10 @@ def map_int_g(self, expr): # {{{ generator - base_kernel = kernel.get_base_kernel() - from sumpy.p2p import P2PMatrixGenerator mat_gen = P2PMatrixGenerator(actx.context, - source_kernels=(base_kernel,), - target_kernels=(target_base_kernel,), + source_kernels=(kernel,), + target_kernels=(expr.target_kernel,), exclude_self=self.exclude_self) # }}} @@ -490,7 +487,7 @@ def map_int_g(self, expr): # {{{ kernel args # NOTE: copied from pytential.symbolic.primitives.IntG - kernel_args = base_kernel.get_args() + base_kernel.get_source_args() + kernel_args = kernel.get_args() + kernel.get_source_args() kernel_args = {arg.loopy_arg.name for arg in kernel_args} kernel_args = _get_layer_potential_args( @@ -638,7 +635,6 @@ def map_int_g(self, expr): expr.target.geometry, expr.target.discr_stage) actx = self.array_context - target_base_kernel = expr.target_kernel.get_base_kernel() result = 0 for kernel, density in zip(expr.source_kernels, expr.densities): @@ -659,12 +655,10 @@ def map_int_g(self, expr): # {{{ generator - base_kernel = kernel.get_base_kernel() - from sumpy.p2p import P2PMatrixSubsetGenerator mat_gen = P2PMatrixSubsetGenerator(actx.context, - source_kernels=(base_kernel,), - target_kernels=(target_base_kernel,), + source_kernels=(kernel,), + target_kernels=(expr.target_kernel,), exclude_self=self.exclude_self) # }}} @@ -674,7 +668,7 @@ def map_int_g(self, expr): # {{{ kernel args # NOTE: copied from pytential.symbolic.primitives.IntG - kernel_args = base_kernel.get_args() + base_kernel.get_source_args() + kernel_args = kernel.get_args() + kernel.get_source_args() kernel_args = {arg.loopy_arg.name for arg in kernel_args} kernel_args = _get_layer_potential_args(