From d9127a3309db56cfdba4bdcb30ed4670c5db521f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 18:23:04 -0700 Subject: [PATCH 01/20] Minor cleanups --- pytential/qbx/__init__.py | 5 ++++- pytential/qbx/geometry.py | 2 +- pytential/qbx/refinement.py | 9 +++++++++ pytential/qbx/utils.py | 8 +++++--- 4 files changed, 19 insertions(+), 5 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 9c9b4ea3..85f53984 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -367,7 +367,10 @@ class QBXLayerPotentialSource(LayerPotentialSource): return self._centers_of_mass_for_discr(self.base_fine_density_discr) def _panel_sizes_for_discr(self, discr, last_dim_length): - assert last_dim_length in ("nsources", "ncenters", "npanels") + if last_dim_length not in ("nsources", "ncenters", "npanels"): + raise ValueError( + "invalid value of last_dim_length: %s" % last_dim_length) + # To get the panel size this does the equivalent of (∫ 1 ds)**(1/dim). # FIXME: Kernel optimizations diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index ce93e851..211e25f9 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -679,7 +679,7 @@ class QBXFMMGeometryData(object): @memoize_method def tree(self): - """Build and return a :class:`boxtree.tree.TreeWithLinkedPointSources` + """Build and return a :class:`boxtree.tree.Tree` for this source with these targets. |cached| diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 9ae54bb6..94d18fb4 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -270,6 +270,8 @@ class RefinerWrangler(object): self.code_container = code_container self.queue = queue + # {{{ check subroutines for conditions 1-3 + def check_expansion_disks_undisturbed_by_sources(self, lpot_source, tree, peer_lists, refine_flags, debug, wait_for=None): @@ -416,6 +418,8 @@ class RefinerWrangler(object): return (out["refine_flags_updated"].get() == 1).all() + # }}} + def build_tree(self, lpot_source, use_base_fine_discr=False): tb = self.code_container.tree_builder() from pytential.qbx.utils import build_tree_with_qbx_metadata @@ -470,6 +474,8 @@ def make_empty_refine_flags(queue, lpot_source, use_base_fine_discr=False): return result +# {{{ main entry point + def refine_for_global_qbx(lpot_source, code_container, group_factory, fine_group_factory, kernel_length_scale=None, # FIXME: Set debug=False once everything works. @@ -621,6 +627,7 @@ def refine_for_global_qbx(lpot_source, code_container, lpot_source = lpot_source.copy(debug=debug, refined_for_global_qbx=True) if len(connections) == 0: + # FIXME: This is inefficient connection = make_same_mesh_connection( lpot_source.density_discr, lpot_source.density_discr) @@ -629,5 +636,7 @@ def refine_for_global_qbx(lpot_source, code_container, return lpot_source, connection +# }}} + # vim: foldmethod=marker:filetype=pyopencl diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index 8b363250..e22336a2 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -150,8 +150,7 @@ def plot_discr(lpot_source, outfilename="discr.pdf"): # }}} -# {{{ tree creation - +# {{{ tree-with-metadata: data structure class TreeWithQBXMetadata(Tree): """A subclass of :class:`boxtree.tree.Tree`. Has all of that class's @@ -229,9 +228,12 @@ class TreeWithQBXMetadata(Tree): """ pass +# }}} -MAX_REFINE_WEIGHT = 64 +# {{{ tree-with-metadata: creation + +MAX_REFINE_WEIGHT = 64 def build_tree_with_qbx_metadata( queue, tree_builder, lpot_source, targets_list=(), -- GitLab From 1357c2db70c1b08046e60478cd93b6d1887d4600 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 18:50:11 -0700 Subject: [PATCH 02/20] Kill some dead code --- pytential/qbx/geometry.py | 61 --------------------------------------- 1 file changed, 61 deletions(-) diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 211e25f9..1c32498f 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -118,35 +118,6 @@ class QBXFMMGeometryCodeGetter(object): self.ambient_dim = ambient_dim self.debug = debug - @property - @memoize_method - def pick_expansion_centers(self): - knl = lp.make_kernel( - """{[dim,k,i]: - 0<=dim Date: Fri, 26 May 2017 19:25:18 -0700 Subject: [PATCH 03/20] Move center getter somewhere more internal --- pytential/qbx/__init__.py | 12 ------------ pytential/qbx/utils.py | 20 ++++++++++++++++++-- pytential/symbolic/matrix.py | 4 +++- test/test_global_qbx.py | 12 ++++++++---- test/test_layer_pot.py | 5 ++++- 5 files changed, 33 insertions(+), 20 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 85f53984..9d851a92 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -435,18 +435,6 @@ class QBXLayerPotentialSource(LayerPotentialSource): return self._panel_sizes_for_discr( self.base_fine_density_discr, last_dim_length) - @memoize_method - def centers(self, sign): - adim = self.density_discr.ambient_dim - dim = self.density_discr.dim - - from pytential import sym, bind - with cl.CommandQueue(self.cl_context) as queue: - nodes = bind(self.density_discr, sym.nodes(adim))(queue) - normals = bind(self.density_discr, sym.normal(adim, dim=dim))(queue) - panel_sizes = self.panel_sizes("nsources").with_queue(queue) - return (nodes + normals * sign * panel_sizes / 2).as_vector(np.object) - @memoize_method def weights_and_area_elements(self): import pytential.symbolic.primitives as p diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index e22336a2..cc1a0049 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -67,6 +67,22 @@ QBX_TREE_MAKO_DEFS = r"""//CL:mako// # }}} +# {{{ compute center array + +def get_centers_on_side(lpot_src, sign): + adim = lpot_src.density_discr.ambient_dim + dim = lpot_src.density_discr.dim + + from pytential import sym, bind + with cl.CommandQueue(lpot_src.cl_context) as queue: + nodes = bind(lpot_src.density_discr, sym.nodes(adim))(queue) + normals = bind(lpot_src.density_discr, sym.normal(adim, dim=dim))(queue) + panel_sizes = lpot_src.panel_sizes("nsources").with_queue(queue) + return (nodes + normals * sign * panel_sizes / 2).as_vector(np.object) + +# }}} + + # {{{ interleaver kernel @memoize @@ -99,8 +115,8 @@ def get_interleaved_centers(queue, lpot_source): next to corresponding exterior centers. """ knl = get_interleaver_kernel(lpot_source.density_discr.real_dtype) - int_centers = lpot_source.centers(-1) - ext_centers = lpot_source.centers(+1) + int_centers = get_centers_on_side(lpot_source, -1) + ext_centers = get_centers_on_side(lpot_source, +1) result = [] wait_for = [] diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 9b4622fd..2012b992 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -183,11 +183,13 @@ class MatrixBuilder(EvaluationMapperBase): assert target_discr is source.density_discr + from pytential.qbx.utils import get_centers_on_side + assert abs(expr.qbx_forced_limit) > 0 _, (mat,) = mat_gen(self.queue, target_discr.nodes(), source.fine_density_discr.nodes(), - source.centers(expr.qbx_forced_limit), + get_centers_on_side(source, expr.qbx_forced_limit), **kernel_args) mat = mat.get() diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index 4013b8f3..d5ca949c 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -107,11 +107,13 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): lpot_source, RefinerCodeContainer(cl_ctx), factory, fine_factory, **refiner_extra_kwargs) + from pytential.qbx.utils import get_centers_on_side + discr_nodes = lpot_source.density_discr.nodes().get(queue) fine_discr_nodes = lpot_source.fine_density_discr.nodes().get(queue) - int_centers = lpot_source.centers(-1) + int_centers = get_centers_on_side(lpot_source, -1) int_centers = np.array([axis.get(queue) for axis in int_centers]) - ext_centers = lpot_source.centers(+1) + ext_centers = get_centers_on_side(lpot_source, +1) ext_centers = np.array([axis.get(queue) for axis in ext_centers]) panel_sizes = lpot_source.panel_sizes("npanels").get(queue) fine_panel_sizes = lpot_source.fine_panel_sizes("npanels").get(queue) @@ -232,8 +234,10 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements): lpot_source, conn = QBXLayerPotentialSource(discr, order).with_refinement() del discr - int_centers = lpot_source.centers(-1) - ext_centers = lpot_source.centers(+1) + from pytential.qbx.utils import get_centers_on_side + + int_centers = get_centers_on_side(lpot_source, -1) + ext_centers = get_centers_on_side(lpot_source, +1) # }}} diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 870cc6f0..af43bb94 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -160,7 +160,10 @@ def test_ellipse_eigenvalues(ctx_getter, ellipse_aspect, mode_nr, qbx_order): if 0: # plot geometry, centers, normals - centers = qbx.centers(density_discr, 1) + + from pytential.qbx.utils import get_centers_on_side + centers = get_centers_on_side(qbx, 1) + nodes_h = nodes.get() centers_h = [centers[0].get(), centers[1].get()] pt.plot(nodes_h[0], nodes_h[1], "x-") -- GitLab From bf90a08cbd011e4f26ae3f8db0bceb1c75329fc2 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 19:27:20 -0700 Subject: [PATCH 04/20] Kill unused find_element_centers --- pytential/qbx/geometry.py | 18 ------------------ 1 file changed, 18 deletions(-) diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 1c32498f..c36664a3 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -118,24 +118,6 @@ class QBXFMMGeometryCodeGetter(object): self.ambient_dim = ambient_dim self.debug = debug - @property - @memoize_method - def find_element_centers(self): - knl = lp.make_kernel( - """{[dim,k,i]: - 0<=dim Date: Fri, 26 May 2017 19:45:09 -0700 Subject: [PATCH 05/20] Kill unused _JumpTermArgumentProvider --- pytential/qbx/__init__.py | 39 --------------------------------------- 1 file changed, 39 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 9d851a92..f0ae9955 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -45,45 +45,6 @@ __doc__ = """ """ -# {{{ jump term interface helper - -class _JumpTermArgumentProvider(object): - def __init__(self, discr, density, ds_direction, side=None): - self.discr = discr - self.density = density - self.ds_direction = ds_direction - self.side = side - - @property - def normal(self): - return self.discr.curve.normals.reshape(2, -1).T - - @property - def tangent(self): - return self.discr.curve.tangents.reshape(2, -1).T - - @property - def src_derivative_dir(self): - return self.ds_direction - - @property - def mean_curvature(self): - return self.discr.curve.curvature.reshape(-1) - - @property - def density_0(self): - return self.density.reshape(-1) - - @property - @memoize_method - def density_0_prime(self): - diff_mat = self.discr.curve.expansion.get_differentiation_matrix() - return (2 * np.dot(diff_mat, self.density.T).T.reshape(-1) - / self.discr.curve.speed.reshape(-1)) - -# }}} - - class LayerPotentialSource(PotentialSource): pass -- GitLab From bd162ae6ded12e6b0d1d874bc9969328af6b0557 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 20:06:40 -0700 Subject: [PATCH 06/20] Seed the Muller solve to eliminate flakiness (Fixes #26 on Gitlab) --- pytential/muller.py | 4 ++-- test/test_muller.py | 6 +++++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/pytential/muller.py b/pytential/muller.py index 22fe5458..5fc6a823 100644 --- a/pytential/muller.py +++ b/pytential/muller.py @@ -25,7 +25,7 @@ THE SOFTWARE. import numpy as np -def muller_deflate(f, n, maxiter=100, eps=1e-14): +def muller_deflate(f, n, maxiter=100, eps=1e-14, z_start=None): """ :arg n: number of zeros sought :returns: (roots, niter) - roots of the given function; number of @@ -52,7 +52,7 @@ def muller_deflate(f, n, maxiter=100, eps=1e-14): niter.append(niter0) while (np.isnan(roots[i]) or niter[i] == maxiter) and miter < 50: - roots0, niter0 = muller(f_deflated, maxiter, eps) + roots0, niter0 = muller(f_deflated, maxiter, eps, z_start=z_start) roots[i] = roots0 niter[i] = niter0 miter = miter+1 diff --git a/test/test_muller.py b/test/test_muller.py index 30d1a47d..fb23e911 100644 --- a/test/test_muller.py +++ b/test/test_muller.py @@ -37,10 +37,14 @@ def test_muller(true_roots): :arg n: number of zeros sought :return: (roots, niter) """ + np.random.seed(15) + z_start = np.random.rand(3) + 1j*np.random.rand(3) + eps = 1e-12 from pytential.muller import muller_deflate roots, niter = muller_deflate( - lambda z: poly_with_roots(z, true_roots), len(true_roots)) + lambda z: poly_with_roots(z, true_roots), len(true_roots), + z_start=z_start) for r_i in roots: min_dist, true_root = min( -- GitLab From f2421e962a78b2fda44e0960ae7cf0e023022e58 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 22:04:13 -0700 Subject: [PATCH 07/20] Untangle exp radius/tunnel radius/panel size, shrink QBXLayerPotentialSource class def --- pytential/qbx/__init__.py | 199 ++++++++++++---------------------- pytential/qbx/geometry.py | 13 ++- pytential/qbx/refinement.py | 10 +- pytential/qbx/target_assoc.py | 49 +++++---- pytential/qbx/utils.py | 128 +++++++++++++++++++++- test/test_global_qbx.py | 22 ++-- 6 files changed, 247 insertions(+), 174 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index f0ae9955..8918e5a2 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -25,7 +25,6 @@ THE SOFTWARE. import six -import loopy as lp import numpy as np from pytools import memoize_method from meshmode.discretization import Discretization @@ -89,13 +88,14 @@ class QBXLayerPotentialSource(LayerPotentialSource): .. attribute :: qbx_order .. attribute :: fmm_order .. attribute :: cl_context - .. automethod :: centers - .. automethod :: panel_sizes .. automethod :: weights_and_area_elements .. automethod :: with_refinement See :ref:`qbxguts` for some information on the inner workings of this. """ + + # {{{ constructor / copy + def __init__(self, density_discr, fine_order, qbx_order=None, fmm_order=None, @@ -179,10 +179,13 @@ class QBXLayerPotentialSource(LayerPotentialSource): else self.refined_for_global_qbx), performance_data_file=self.performance_data_file) + # }}} + @property def base_fine_density_discr(self): - """The fine density discretization before upsampling is applied. + """The refined, interpolation-focused density discretization (no oversampling). """ + # FIXME: Maybe rename interp_refined_discr return (self._base_resampler.to_discr if self._base_resampler is not None else self.density_discr) @@ -190,6 +193,9 @@ class QBXLayerPotentialSource(LayerPotentialSource): @property @memoize_method def fine_density_discr(self): + """The refined, quadrature-focused density discretization (with upsampling). + """ + # FIXME: Maybe rename quad_refined_discr from meshmode.discretization.poly_element import ( QuadratureSimplexGroupFactory) @@ -212,22 +218,6 @@ class QBXLayerPotentialSource(LayerPotentialSource): return conn - def el_view(self, discr, group_nr, global_array): - """Return a view of *global_array* of shape - ``(..., discr.groups[group_nr].nelements)`` - where *global_array* is of shape ``(..., nelements)``, - where *nelements* is the global (per-discretization) node count. - """ - - group = discr.groups[group_nr] - el_nr_base = sum(group.nelements for group in discr.groups[:group_nr]) - - return global_array[ - ..., el_nr_base:el_nr_base + group.nelements] \ - .reshape( - global_array.shape[:-1] - + (group.nelements,)) - @property @memoize_method def refiner_code_container(self): @@ -266,7 +256,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): @memoize_method def h_max(self): with cl.CommandQueue(self.cl_context) as queue: - panel_sizes = self.panel_sizes("npanels").with_queue(queue) + panel_sizes = self._panel_sizes("npanels").with_queue(queue) return np.asscalar(cl.array.max(panel_sizes).get()) @property @@ -289,112 +279,72 @@ class QBXLayerPotentialSource(LayerPotentialSource): def complex_dtype(self): return self.density_discr.complex_dtype - def _centers_of_mass_for_discr(self, discr): - knl = lp.make_kernel( - """{[dim,k,i]: - 0<=dim= h / 2, \ - (dist, h, centers_panel.element_nr, sources_panel.element_nr) + rad = expansion_radii[centers_panel.element_nr] + assert dist >= rad, \ + (dist, rad, centers_panel.element_nr, sources_panel.element_nr) def check_sufficient_quadrature_resolution(centers_panel, sources_panel): h = fine_panel_sizes[sources_panel.element_nr] @@ -247,7 +246,8 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements): rng = PhiloxGenerator(cl_ctx, seed=RNG_SEED) nsources = lpot_source.density_discr.nnodes noise = rng.uniform(queue, nsources, dtype=np.float, a=0.01, b=1.0) - panel_sizes = lpot_source.panel_sizes("nsources").with_queue(queue) + tunnel_radius = \ + lpot_source._close_target_tunnel_radius("nsources").with_queue(queue) def targets_from_sources(sign, dist): from pytential import sym, bind @@ -258,8 +258,8 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements): from pytential.target import PointsTarget - int_targets = PointsTarget(targets_from_sources(-1, noise * panel_sizes / 2)) - ext_targets = PointsTarget(targets_from_sources(+1, noise * panel_sizes / 2)) + int_targets = PointsTarget(targets_from_sources(-1, noise * tunnel_radius)) + ext_targets = PointsTarget(targets_from_sources(+1, noise * tunnel_radius)) far_targets = PointsTarget(targets_from_sources(+1, FAR_TARGET_DIST_FROM_SOURCE)) # Create target discretizations. @@ -294,7 +294,7 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements): QBXTargetAssociator(cl_ctx)(lpot_source, target_discrs) .get(queue=queue)) - panel_sizes = lpot_source.panel_sizes("nsources").get(queue) + expansion_radii = lpot_source._expansion_radii("nsources").get(queue) int_centers = np.array([axis.get(queue) for axis in int_centers]) ext_centers = np.array([axis.get(queue) for axis in ext_centers]) @@ -314,7 +314,7 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements): target_to_source_result, target_to_side_result): assert (target_to_side_result == true_side).all() dists = la.norm((targets.T - centers.T[target_to_source_result]), axis=1) - assert (dists <= panel_sizes[target_to_source_result] / 2).all() + assert (dists <= expansion_radii[target_to_source_result]).all() # Checks that far targets are not assigned a center. def check_far_targets(target_to_source_result): -- GitLab From 11ee435351e40e705325cc6df59b40413ff71413 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 22:13:35 -0700 Subject: [PATCH 08/20] Re-add missing panel_sizes assignment to refinement test --- test/test_global_qbx.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index 1e394595..7326efc1 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -116,6 +116,7 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): ext_centers = get_centers_on_side(lpot_source, +1) ext_centers = np.array([axis.get(queue) for axis in ext_centers]) expansion_radii = lpot_source._expansion_radii("npanels").get(queue) + panel_sizes = lpot_source._panel_sizes("npanels").get(queue) fine_panel_sizes = lpot_source._fine_panel_sizes("npanels").get(queue) # {{{ check if satisfying criteria -- GitLab From 792c943a6980e9ced05fbe443666d86e8deb3466 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 22:18:01 -0700 Subject: [PATCH 09/20] Flake8 whitespace fix --- pytential/qbx/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 8918e5a2..70512a10 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -365,7 +365,6 @@ class QBXLayerPotentialSource(LayerPotentialSource): return (area_element.with_queue(queue)*qweight).with_queue(None) - # {{{ helpers for symbolic operator processing def preprocess_optemplate(self, name, discretizations, expr): -- GitLab From 5d17d57d17332c4edd78bb09c8b753afe4c0da63 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 22:32:49 -0700 Subject: [PATCH 10/20] Kill more dead code --- pytential/qbx/geometry.py | 26 -------------------------- 1 file changed, 26 deletions(-) diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 915f55b6..cb3ac634 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -461,7 +461,6 @@ class QBXFMMGeometryData(object): .. automethod:: global_qbx_centers() .. automethod:: user_target_to_center() .. automethod:: center_to_tree_targets() - .. automethod:: global_qbx_centers_box_target_lists() .. automethod:: non_qbx_box_target_lists() .. automethod:: plot() """ @@ -849,31 +848,6 @@ class QBXFMMGeometryData(object): return result - @memoize_method - def global_qbx_centers_box_target_lists(self): - """Build a list of targets per box consisting only of global QBX centers. - Returns a :class:`boxtree.tree.FilteredTargetListsInUserOrder`. - (I.e. no new target order is created for these targets, as we expect - there to be (relatively) not many of them.) - - |cached| - """ - - center_info = self.center_info() - with cl.CommandQueue(self.cl_context) as queue: - logger.info("find global qbx centers box target list: start") - - flags = cl.array.zeros(queue, self.tree().ntargets, np.int8) - - flags[:center_info.ncenters] = self.global_qbx_flags() - - from boxtree.tree import filter_target_lists_in_user_order - result = filter_target_lists_in_user_order(queue, self.tree(), flags) - - logger.info("find global qbx centers box target list: done") - - return result.with_queue(None) - @memoize_method def non_qbx_box_target_lists(self): """Build a list of targets per box that don't need to bother with QBX. -- GitLab From d1554561de3d13e0b7f41d504d25b0b424320ecd Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 22:36:22 -0700 Subject: [PATCH 11/20] Mop up another orphaned centers computation --- pytential/qbx/__init__.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 70512a10..1cf8122f 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -658,6 +658,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): strengths = (evaluate(insn.density).with_queue(queue) * self.weights_and_area_elements()) + import pytential.qbx.utils as utils + # FIXME: Do this all at once result = [] for o in insn.outputs: @@ -673,7 +675,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): evt, output_for_each_kernel = lpot_applier( queue, target_discr.nodes(), self.fine_density_discr.nodes(), - self.centers(o.qbx_forced_limit), + utils.get_centers_on_side(self, o.qbx_forced_limit), [strengths], **kernel_args) result.append((o.name, output_for_each_kernel[o.kernel_index])) else: -- GitLab From adaf1770ed72ec0d14d3754f23d84c364b49165d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 23:04:29 -0700 Subject: [PATCH 12/20] Shuffle utils order --- pytential/qbx/utils.py | 46 +++++++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index 702fda66..27436461 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -181,6 +181,29 @@ def get_centers_on_side(lpot_src, sign): # }}} +# {{{ el_view + +def el_view(discr, group_nr, global_array): + """Return a view of *global_array* of shape + ``(..., discr.groups[group_nr].nelements)`` + where *global_array* is of shape ``(..., nelements)``, + where *nelements* is the global (per-discretization) node count. + """ + + group = discr.groups[group_nr] + el_nr_base = sum(group.nelements for group in discr.groups[:group_nr]) + + return global_array[ + ..., el_nr_base:el_nr_base + group.nelements] \ + .reshape( + global_array.shape[:-1] + + (group.nelements,)) + +# }}} + + +# {{{ make interleaved centers + # {{{ interleaver kernel @memoize @@ -205,29 +228,6 @@ def get_interleaver_kernel(dtype): # }}} -# {{{ el_view - -def el_view(discr, group_nr, global_array): - """Return a view of *global_array* of shape - ``(..., discr.groups[group_nr].nelements)`` - where *global_array* is of shape ``(..., nelements)``, - where *nelements* is the global (per-discretization) node count. - """ - - group = discr.groups[group_nr] - el_nr_base = sum(group.nelements for group in discr.groups[:group_nr]) - - return global_array[ - ..., el_nr_base:el_nr_base + group.nelements] \ - .reshape( - global_array.shape[:-1] - + (group.nelements,)) - -# }}} - - -# {{{ make interleaved centers - def get_interleaved_centers(queue, lpot_source): """ Return an array of shape (dim, ncenters) in which interior centers are placed -- GitLab From 4941d09eab95d480c146f8717fb6059a0dd13c3e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 23:11:54 -0700 Subject: [PATCH 13/20] Get rid of CenterInfo --- pytential/qbx/__init__.py | 9 ++- pytential/qbx/fmm.py | 16 ++--- pytential/qbx/geometry.py | 143 +++++++++++++------------------------- 3 files changed, 61 insertions(+), 107 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 1cf8122f..8fcd8c97 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -540,7 +540,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): # }}} - if len(geo_data.global_qbx_centers()) != geo_data.center_info().ncenters: + if len(geo_data.global_qbx_centers()) != geo_data.ncenters: raise NotImplementedError("geometry has centers requiring local QBX") from pytential.qbx.geometry import target_state @@ -699,12 +699,11 @@ class QBXLayerPotentialSource(LayerPotentialSource): (target_discr, qbx_forced_limit), )) - # center_info is independent of targets - center_info = geo_data.center_info() + # center-related info is independent of targets # First ncenters targets are the centers tgt_to_qbx_center = ( - geo_data.user_target_to_center()[center_info.ncenters:] + geo_data.user_target_to_center()[geo_data.ncenters:] .copy(queue=queue)) qbx_tgt_numberer = self.get_qbx_target_numberer( @@ -737,7 +736,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): queue, targets=target_discr.nodes(), sources=self.fine_density_discr.nodes(), - centers=center_info.centers, + centers=geo_data.centers, strengths=[strengths], qbx_tgt_numbers=qbx_tgt_numbers, qbx_center_numbers=qbx_center_numbers, diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index e3508b70..7fd7f49d 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -153,7 +153,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), return cl.array.zeros( self.queue, - (self.geo_data.center_info().ncenters, + (self.geo_data.ncenters, len(qbx_l_expn)), dtype=self.dtype) @@ -212,7 +212,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), self.queue, global_qbx_centers=geo_data.global_qbx_centers(), qbx_center_to_target_box=geo_data.qbx_center_to_target_box(), - qbx_centers=geo_data.center_info().centers, + qbx_centers=geo_data.centers(), source_box_starts=starts, source_box_lists=lists, @@ -230,7 +230,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), qbx_expansions = self.qbx_local_expansion_zeros() geo_data = self.geo_data - if geo_data.center_info().ncenters == 0: + if geo_data.ncenters == 0: return qbx_expansions traversal = geo_data.traversal() @@ -249,7 +249,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), qbx_center_to_target_box=geo_data.qbx_center_to_target_box(), centers=self.tree.box_centers, - qbx_centers=geo_data.center_info().centers, + qbx_centers=geo_data.centers(), src_expansions=source_mpoles_view, src_base_ibox=source_level_start_ibox, @@ -273,7 +273,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), qbx_expansions = self.qbx_local_expansion_zeros() geo_data = self.geo_data - if geo_data.center_info().ncenters == 0: + if geo_data.ncenters == 0: return qbx_expansions trav = geo_data.traversal() @@ -294,7 +294,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), target_base_ibox=target_level_start_ibox, centers=self.tree.box_centers, - qbx_centers=geo_data.center_info().centers, + qbx_centers=geo_data.centers(), expansions=target_locals_view, qbx_expansions=qbx_expansions, @@ -322,7 +322,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), qbxl2p = self.code.qbxl2p(self.qbx_order) evt, pot_res = qbxl2p(self.queue, - qbx_centers=geo_data.center_info().centers, + qbx_centers=geo_data.centers(), global_qbx_centers=geo_data.global_qbx_centers(), center_to_targets_starts=ctt.starts, @@ -652,7 +652,7 @@ def write_performance_model(outf, geo_data): qbx_center_to_target_box = geo_data.qbx_center_to_target_box() center_to_targets_starts = geo_data.center_to_tree_targets().starts - ncenters = geo_data.center_info().ncenters + ncenters = geo_data.ncenters outf.write("ncenters = {cost}\n" .format(cost=ncenters)) diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index cb3ac634..2f29bb7b 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -353,32 +353,6 @@ class TargetInfo(DeviceDataRecord): """ -class CenterInfo(DeviceDataRecord): - """Information on location of QBX centers. - Returned from :meth:`QBXFMMGeometryData.center_info`. - - .. attribute:: centers - - Shape: ``[dim][ncenters]`` - - .. attribute:: sides - - Shape: ``[ncenters]`` - - -1 for inside, +1 for outside (relative to normal) - - .. attribute:: radii - - Shape: ``[ncenters]`` - - .. attribute:: ncenters - """ - - @property - def ncenters(self): - return len(self.centers[0]) - - class CenterToTargetList(DeviceDataRecord): """A lookup table of targets covered by each QBX disk. Indexed by global number of QBX center, ``lists[start[i]:start[i+1]]`` indicates numbers @@ -452,7 +426,9 @@ class QBXFMMGeometryData(object): .. rubric :: Methods - .. automethod:: center_info() + .. attribute:: ncenters + + .. automethod:: centers() .. automethod:: target_info() .. automethod:: tree() .. automethod:: traversal() @@ -493,19 +469,10 @@ class QBXFMMGeometryData(object): def coord_dtype(self): return self.lpot_source.fine_density_discr.nodes().dtype - # {{{ center info - - @memoize_method - def kept_center_indices(self, el_group): - """Return indices of unit nodes (of the target discretization) - that will give rise to QBX centers. - """ - - # FIXME: Be more careful about which nodes to keep - return np.arange(0, el_group.nunit_nodes) + # {{{ centers @memoize_method - def center_info(self): + def centers(self): """ Return a :class:`CenterInfo`. |cached| """ @@ -513,16 +480,14 @@ class QBXFMMGeometryData(object): with cl.CommandQueue(self.cl_context) as queue: from pytential.qbx.utils import get_interleaved_centers - centers = get_interleaved_centers(queue, lpot_source) - # sides = cl.array.arange(queue, len(centers[0]), dtype=np.int32) - # sides = 2 * (sides & 1) - 1 - # radii = lpot_source.panel_sizes("ncenters").with_queue(queue) / 2 + from pytools.obj_array import make_obj_array + return make_obj_array([ + ccomp.with_queue(None) + for ccomp in get_interleaved_centers(queue, lpot_source)]) - return CenterInfo( - # FIXME: Conceivably unused - # sides=sides, - # radii=radii, - centers=centers).with_queue(None) + @property + def ncenters(self): + return len(self.centers()[0]) # }}} @@ -536,9 +501,7 @@ class QBXFMMGeometryData(object): lpot_src = self.lpot_source with cl.CommandQueue(self.cl_context) as queue: - center_info = self.center_info() - - ntargets = center_info.ncenters + ntargets = self.ncenters target_discr_starts = [] for target_discr, qbx_side in self.target_discrs_and_qbx_sides: @@ -552,8 +515,8 @@ class QBXFMMGeometryData(object): self.coord_dtype) code_getter.copy_targets_kernel()( queue, - targets=targets[:, :center_info.ncenters], - points=center_info.centers) + targets=targets[:, :self.ncenters], + points=self.centers()) for start, (target_discr, _) in zip( target_discr_starts, self.target_discrs_and_qbx_sides): @@ -567,20 +530,19 @@ class QBXFMMGeometryData(object): target_discr_starts=target_discr_starts, ntargets=ntargets).with_queue(None) - def target_radii(self): + # FIXME unused + def xxx_target_radii(self): """Shape: ``[ntargets]`` A list of target radii for FMM tree construction. In this list, the QBX - centers have the radii determined by :meth:`center_info`, and all other - targets have radius zero. + centers have nonzero radii, and all other targets have radius zero. """ tgt_info = self.target_info() - center_info = self.center_info() with cl.CommandQueue(self.cl_context) as queue: target_radii = cl.array.zeros(queue, tgt_info.ntargets, self.coord_dtype) - target_radii[:center_info.ncenters] = center_info.radii + # target_radii[:self.ncenters] = center_info.radii return target_radii.with_queue(None) @@ -591,12 +553,11 @@ class QBXFMMGeometryData(object): Shape: ``[ntargets]``, dtype: int8""" tgt_info = self.target_info() - center_info = self.center_info() with cl.CommandQueue(self.cl_context) as queue: target_side_preferences = cl.array.empty( queue, tgt_info.ntargets, np.int8) - target_side_preferences[:center_info.ncenters] = 0 + target_side_preferences[:self.ncenters] = 0 for tdstart, (target_discr, qbx_side) in \ zip(tgt_info.target_discr_starts, @@ -668,15 +629,14 @@ class QBXFMMGeometryData(object): @memoize_method def qbx_center_to_target_box(self): - """Return a lookup table of length :attr:`CenterInfo.ncenters` - (see :meth:`center_info`) indicating the target box in which each + """Return a lookup table of length :attr:`ncenters` + indicating the target box in which each QBX disk is located. |cached| """ tree = self.tree() trav = self.traversal() - center_info = self.center_info() qbx_center_to_target_box_lookup = \ self.code_getter.qbx_center_to_target_box_lookup( @@ -709,18 +669,17 @@ class QBXFMMGeometryData(object): box_target_starts=tree.box_target_starts, box_target_counts_nonchild=tree.box_target_counts_nonchild, user_target_from_tree_target=user_target_from_tree_target, - ncenters=center_info.ncenters) + ncenters=self.ncenters) return qbx_center_to_target_box.with_queue(None) @memoize_method def global_qbx_flags(self): """Return an array of :class:`numpy.int8` of length - :attr:`CenterInfo.ncenters` (see :meth:`center_info`) indicating - whether each center can use gloal QBX, i.e. whether a single expansion - can mediate interactions from *all* sources to all targets for which it - is valid. If global QBX can be used, the center's entry will be 1, - otherwise it will be 0. + :attr:`ncenters` indicating whether each center can use gloal QBX, i.e. + whether a single expansion can mediate interactions from *all* sources + to all targets for which it is valid. If global QBX can be used, the + center's entry will be 1, otherwise it will be 0. (If not, local QBX is needed, and the center may only be able to mediate some of the interactions to a given target.) @@ -728,10 +687,8 @@ class QBXFMMGeometryData(object): |cached| """ - center_info = self.center_info() - with cl.CommandQueue(self.cl_context) as queue: - result = cl.array.empty(queue, center_info.ncenters, np.int8) + result = cl.array.empty(queue, self.ncenters, np.int8) result.fill(1) return result.with_queue(None) @@ -750,7 +707,7 @@ class QBXFMMGeometryData(object): logger.info("find global qbx centers: start") result, count, _ = copy_if( - cl.array.arange(queue, self.center_info().ncenters, + cl.array.arange(queue, self.ncenters, tree.particle_id_dtype), "global_qbx_flags[i] != 0", extra_args=[ @@ -778,16 +735,15 @@ class QBXFMMGeometryData(object): tgt_assoc = QBXTargetAssociator(self.cl_context) tgt_info = self.target_info() - center_info = self.center_info() from pytential.target import PointsTarget with cl.CommandQueue(self.cl_context) as queue: target_side_prefs = (self - .target_side_preferences()[center_info.ncenters:].get(queue=queue)) + .target_side_preferences()[self.ncenters:].get(queue=queue)) target_discrs_and_qbx_sides = [( - PointsTarget(tgt_info.targets[:, center_info.ncenters:]), + PointsTarget(tgt_info.targets[:, self.ncenters:]), target_side_prefs.astype(np.int32))] # FIXME: try block... @@ -799,8 +755,8 @@ class QBXFMMGeometryData(object): with cl.CommandQueue(self.cl_context) as queue: result = cl.array.empty(queue, tgt_info.ntargets, tree.particle_id_dtype) - result[:center_info.ncenters].fill(target_state.NO_QBX_NEEDED) - result[center_info.ncenters:] = tgt_assoc_result.target_to_center + result[:self.ncenters].fill(target_state.NO_QBX_NEEDED) + result[self.ncenters:] = tgt_assoc_result.target_to_center return result @@ -812,7 +768,6 @@ class QBXFMMGeometryData(object): |cached| """ - center_info = self.center_info() user_ttc = self.user_target_to_center() with cl.CommandQueue(self.cl_context) as queue: @@ -838,7 +793,7 @@ class QBXFMMGeometryData(object): center_target_starts, targets_sorted_by_center, _ = \ self.code_getter.key_value_sort(queue, filtered_tree_ttc, filtered_target_ids, - center_info.ncenters, tree_ttc.dtype) + self.ncenters, tree_ttc.dtype) logger.info("build center -> targets lookup table: done") @@ -869,7 +824,7 @@ class QBXFMMGeometryData(object): # 'flags' is in user order, and should be. - nqbx_centers = self.center_info().ncenters + nqbx_centers = self.ncenters flags[:nqbx_centers] = 0 from boxtree.tree import filter_target_lists_in_tree_order @@ -908,22 +863,22 @@ class QBXFMMGeometryData(object): # {{{ draw centers and circles - center_info = self.center_info() + centers = self.centers() centers = [ - center_info.centers[0].get(queue), - center_info.centers[1].get(queue)] + centers[0].get(queue), + centers[1].get(queue)] pt.plot(centers[0][global_flags == 0], centers[1][global_flags == 0], "oc", label="centers needing local qbx") ax = pt.gca() - for icenter, (cx, cy, r) in enumerate(zip( - centers[0], centers[1], center_info.radii.get(queue))): - ax.add_artist( - pt.Circle((cx, cy), r, fill=False, ls="dotted", lw=1)) - pt.text(cx, cy, - str(icenter), fontsize=8, - ha="left", va="center", - bbox=dict(facecolor='white', alpha=0.5, lw=0)) + # for icenter, (cx, cy, r) in enumerate(zip( + # centers[0], centers[1], center_info.radii.get(queue))): + # ax.add_artist( + # pt.Circle((cx, cy), r, fill=False, ls="dotted", lw=1)) + # pt.text(cx, cy, + # str(icenter), fontsize=8, + # ha="left", va="center", + # bbox=dict(facecolor='white', alpha=0.5, lw=0)) # }}} @@ -950,9 +905,9 @@ class QBXFMMGeometryData(object): tccount = 0 checked = 0 for tx, ty, tcenter in zip( - targets[0][center_info.ncenters:], - targets[1][center_info.ncenters:], - ttc[center_info.ncenters:]): + targets[0][self.ncenters:], + targets[1][self.ncenters:], + ttc[self.ncenters:]): checked += 1 if tcenter >= 0: tccount += 1 -- GitLab From b312cdd252e9a7f9bdfd0add50e1aa36272933ea Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 23:23:12 -0700 Subject: [PATCH 14/20] Clean up ReST doc markup in TreeWithQBXMetadata --- pytential/qbx/utils.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index 27436461..b29f11ab 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -304,41 +304,41 @@ class TreeWithQBXMetadata(Tree): .. attribute:: box_to_qbx_panel_starts - ``box_id_t [nboxes + 1]`` + ``box_id_t [nboxes + 1]`` .. attribute:: box_to_qbx_panel_lists - ``particle_id_t [*]`` + ``particle_id_t [*]`` .. rubric:: Box to QBX sources .. attribute:: box_to_qbx_source_starts - ``box_id_t [nboxes + 1]`` + ``box_id_t [nboxes + 1]`` .. attribute:: box_to_qbx_source_lists - ``particle_id_t [*]`` + ``particle_id_t [*]`` .. rubric:: Box to QBX centers .. attribute:: box_to_qbx_center_starts - ``box_id_t [nboxes + 1]`` + ``box_id_t [nboxes + 1]`` .. attribute:: box_to_qbx_center_lists - ``particle_id_t [*]`` + ``particle_id_t [*]`` .. rubric:: Box to QBX targets .. attribute:: box_to_qbx_target_starts - ``box_id_t [nboxes + 1]`` + ``box_id_t [nboxes + 1]`` .. attribute:: box_to_qbx_target_lists - ``particle_id_t [*]`` + ``particle_id_t [*]`` .. ------------------------------------------------------------------------ .. rubric:: Panel properties @@ -346,11 +346,11 @@ class TreeWithQBXMetadata(Tree): .. attribute:: qbx_panel_to_source_starts - ``particle_id_t [nqbxpanels + 1]`` + ``particle_id_t [nqbxpanels + 1]`` .. attribute:: qbx_panel_to_center_starts - ``particle_id_t [nqbxpanels + 1]`` + ``particle_id_t [nqbxpanels + 1]`` .. ------------------------------------------------------------------------ .. rubric:: Particle order indices -- GitLab From d02fea4e5353aa286f22ac3929ce99a61feedf8d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 23:23:41 -0700 Subject: [PATCH 15/20] Delete unused qbx_center_for_global_tester --- pytential/qbx/geometry.py | 98 --------------------------------------- 1 file changed, 98 deletions(-) diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 2f29bb7b..3a5bb324 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -197,104 +197,6 @@ class QBXFMMGeometryCodeGetter(object): from boxtree.area_query import LeavesToBallsLookupBuilder return LeavesToBallsLookupBuilder(self.cl_context) - # {{{ check if a center may be used with global QBX - - @memoize_method - def qbx_center_for_global_tester(self, - coord_dtype, box_id_dtype, particle_id_dtype): - from pyopencl.elementwise import ElementwiseTemplate - return ElementwiseTemplate( - arguments="""//CL:mako// - /* input */ - %for iaxis in range(ambient_dim): - coord_t *center_${iaxis}, /* [ncenters] */ - %endfor - coord_t *radii, /* [ncenters] */ - box_id_t *qbx_center_to_target_box, /* [ncenters] */ - - box_id_t *neighbor_source_boxes_starts, - box_id_t *neighbor_source_boxes_lists, - particle_id_t *box_point_source_starts, - particle_id_t *box_point_source_counts_cumul, - %for iaxis in range(ambient_dim): - coord_t *point_sources_${iaxis}, - %endfor - - /* output */ - char *center_may_use_global_qbx, - """, - operation=r"""//CL:mako// - particle_id_t icenter = i; - - %for iaxis in range(ambient_dim): - coord_t my_center_${iaxis} = center_${iaxis}[icenter]; - %endfor - coord_t radius = radii[icenter]; - - // {{{ see if there are sources close enough to require local QBX - - bool found_too_close = false; - - coord_t radius_squared = radius * radius; - - box_id_t itgt_box = qbx_center_to_target_box[icenter]; - - box_id_t nb_src_start = neighbor_source_boxes_starts[itgt_box]; - box_id_t nb_src_stop = neighbor_source_boxes_starts[itgt_box+1]; - - for ( - box_id_t ibox_list = nb_src_start; - ibox_list < nb_src_stop && !found_too_close; - ++ibox_list) - { - box_id_t neighbor_box_id = - neighbor_source_boxes_lists[ibox_list]; - - box_id_t bps_start = box_point_source_starts[neighbor_box_id]; - box_id_t bps_stop = - bps_start + box_point_source_counts_cumul[neighbor_box_id]; - for ( - box_id_t ipsrc = bps_start; - ipsrc < bps_stop; - ++ipsrc) - { - %for iaxis in range(ambient_dim): - coord_t psource_${iaxis} = point_sources_${iaxis}[ipsrc]; - %endfor - - coord_t dist_squared = 0 - %for iaxis in range(ambient_dim): - + (my_center_${iaxis} - psource_${iaxis}) - * (my_center_${iaxis} - psource_${iaxis}) - %endfor - ; - - const coord_t too_close_slack = (1+1e-2)*(1+1e-2); - found_too_close = found_too_close - || (dist_squared*too_close_slack < radius_squared); - - if (found_too_close) - break; - } - } - - // }}} - - center_may_use_global_qbx[icenter] = !found_too_close; - """, - name="qbx_test_center_for_global").build( - self.cl_context, - type_aliases=( - ("box_id_t", box_id_dtype), - ("particle_id_t", particle_id_dtype), - ("coord_t", coord_dtype), - ), - var_values=( - ("ambient_dim", self.ambient_dim), - )) - - # }}} - @property @memoize_method def key_value_sort(self): -- GitLab From ab9871dd2d589b4f5c4d0c4b8895307d42678099 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 27 May 2017 07:00:55 -0700 Subject: [PATCH 16/20] Add missing call parens for centers --- pytential/qbx/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 8fcd8c97..759166d5 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -736,7 +736,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): queue, targets=target_discr.nodes(), sources=self.fine_density_discr.nodes(), - centers=geo_data.centers, + centers=geo_data.centers(), strengths=[strengths], qbx_tgt_numbers=qbx_tgt_numbers, qbx_center_numbers=qbx_center_numbers, -- GitLab From 2e70a94ade934d33c4c891f3eacadd004d947d6b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 27 May 2017 12:55:09 -0700 Subject: [PATCH 17/20] Improve target_assoc comments/panel size quantity usage --- pytential/qbx/target_assoc.py | 38 ++++++++++++++++++++++++++--------- 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/pytential/qbx/target_assoc.py b/pytential/qbx/target_assoc.py index 6bc97ad3..4255d865 100644 --- a/pytential/qbx/target_assoc.py +++ b/pytential/qbx/target_assoc.py @@ -398,16 +398,29 @@ class QBXTargetAssociator(object): # Perform a space invader query over the sources. source_slice = tree.user_source_ids[tree.qbx_user_source_slice] sources = [axis.with_queue(queue)[source_slice] for axis in tree.sources] - panel_sizes = lpot_source._panel_sizes("nsources").with_queue(queue) - - # FIXME: (AK) Not sure what's going on -- need comment + tunnel_radius_by_source = \ + lpot_source._close_target_tunnel_radius("nsources").with_queue(queue) + + # Target-marking algorithm (TGTMARK): + # + # (1) Use a space invader query to tag each leaf box that intersects with the + # "near-source-detection tunnel" with the distance to the closest source. + # + # (2) Do an area query around all targets with the radius resulting + # from the space invader query, enumerate sources in that vicinity. + # If a source is found whose distance to the target is less than the + # source's tunnel radius, mark that target as pending. + # (or below: mark the source for refinement) + + # Note that this comment is referred to below by "TGTMARK". If you + # remove this comment or change the algorithm here, make sure that + # the reference below is still accurate. box_to_search_dist, evt = self.space_invader_query( queue, tree, sources, - # FIXME: Is this the right quantity? (as opposed to expansion radii?) - panel_sizes / 2, + tunnel_radius_by_source, peer_lists, wait_for=wait_for) wait_for = [evt] @@ -473,6 +486,12 @@ class QBXTargetAssociator(object): expansion_radii_by_center_with_stick_out = \ expansion_radii_by_center * (1 + stick_out_factor) + # Idea: + # + # (1) Tag leaf boxes around centers with max distance to usable center. + # (2) Area query from targets with those radii to find closest eligible + # center. + box_to_search_dist, evt = self.space_invader_query( queue, tree, @@ -542,17 +561,16 @@ class QBXTargetAssociator(object): # Perform a space invader query over the sources. source_slice = tree.user_source_ids[tree.qbx_user_source_slice] sources = [axis.with_queue(queue)[source_slice] for axis in tree.sources] - panel_sizes = lpot_source._panel_sizes("nsources").with_queue(queue) + tunnel_radius_by_source = \ + lpot_source._close_target_tunnel_radius("nsources").with_queue(queue) - # FIXME: (AK) Not sure what's going on -- need comment + # See (TGTMARK) above for algorithm. box_to_search_dist, evt = self.space_invader_query( queue, tree, sources, - # FIXME: Is this the right quantity? - # (as opposed to expansion/tunnel radii?) - panel_sizes / 2, + tunnel_radius_by_source, peer_lists, wait_for=wait_for) wait_for = [evt] -- GitLab From b20b624315edde9baec35ee9cb9f520c60e1dbd8 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 27 May 2017 13:21:08 -0700 Subject: [PATCH 18/20] Minor comment add --- pytential/qbx/geometry.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 3a5bb324..b0731d8b 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -215,7 +215,10 @@ class QBXFMMGeometryCodeGetter(object): VectorArg(particle_id_dtype, "filtered_target_id"), VectorArg(particle_id_dtype, "count"), ], + + # "Does this target have a QBX center?" input_expr="(target_to_center[i] >= 0) ? 1 : 0", + scan_expr="a+b", neutral="0", output_statement=""" if (prev_item != item) -- GitLab From d5fc38cb0b05467442ebfff20fd821a18fc1048b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 27 May 2017 16:51:13 -0700 Subject: [PATCH 19/20] Update some renaming FIXMEs --- pytential/qbx/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 759166d5..1fe8d734 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -185,7 +185,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): def base_fine_density_discr(self): """The refined, interpolation-focused density discretization (no oversampling). """ - # FIXME: Maybe rename interp_refined_discr + # FIXME: Maybe rename refined_interp_density_discr return (self._base_resampler.to_discr if self._base_resampler is not None else self.density_discr) @@ -195,7 +195,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): def fine_density_discr(self): """The refined, quadrature-focused density discretization (with upsampling). """ - # FIXME: Maybe rename quad_refined_discr + # FIXME: Maybe rename refined_quad_density_discr from meshmode.discretization.poly_element import ( QuadratureSimplexGroupFactory) -- GitLab From 46d65fe3417eb1ae1d21a5567373b6df6d662a60 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 27 May 2017 16:54:54 -0700 Subject: [PATCH 20/20] Add (optional) flag to enable zero-box-slack FMM --- pytential/qbx/__init__.py | 5 + pytential/qbx/geometry.py | 64 +- pytential/qbx/utils.py | 121 ++- .../symbolic/pde/maxwell/waveguide.py.orig | 989 ++++++++++++++++++ 4 files changed, 1101 insertions(+), 78 deletions(-) create mode 100644 pytential/symbolic/pde/maxwell/waveguide.py.orig diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 1fe8d734..0d058c82 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -107,6 +107,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): # FIXME default debug=False once everything works debug=True, refined_for_global_qbx=False, + expansion_disks_in_tree_have_extent=False, performance_data_file=None): """ :arg fine_order: The total degree to which the (upsampled) @@ -145,6 +146,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): self.debug = debug self.refined_for_global_qbx = refined_for_global_qbx + self.expansion_disks_in_tree_have_extent = \ + expansion_disks_in_tree_have_extent self.performance_data_file = performance_data_file def copy( @@ -177,6 +180,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): refined_for_global_qbx=( refined_for_global_qbx if refined_for_global_qbx is not None else self.refined_for_global_qbx), + expansion_disks_in_tree_have_extent=( + self.expansion_disks_in_tree_have_extent), performance_data_file=self.performance_data_file) # }}} diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index b0731d8b..0eee50d6 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -329,11 +329,14 @@ class QBXFMMGeometryData(object): .. attribute:: coord_dtype - .. rubric :: Methods + .. rubric :: Expansion centers .. attribute:: ncenters - .. automethod:: centers() + .. automethod:: radii() + + .. rubric :: Methods + .. automethod:: target_info() .. automethod:: tree() .. automethod:: traversal() @@ -374,25 +377,36 @@ class QBXFMMGeometryData(object): def coord_dtype(self): return self.lpot_source.fine_density_discr.nodes().dtype - # {{{ centers + # {{{ centers/radii + + @property + def ncenters(self): + return len(self.centers()[0]) @memoize_method def centers(self): - """ Return a :class:`CenterInfo`. |cached| - """ + """ Return an object array of (interleaved) center coordinates. - lpot_source = self.lpot_source + ``coord_t [ambient_dim][ncenters]`` + """ with cl.CommandQueue(self.cl_context) as queue: from pytential.qbx.utils import get_interleaved_centers from pytools.obj_array import make_obj_array return make_obj_array([ ccomp.with_queue(None) - for ccomp in get_interleaved_centers(queue, lpot_source)]) + for ccomp in get_interleaved_centers(queue, self.lpot_source)]) - @property - def ncenters(self): - return len(self.centers()[0]) + @memoize_method + def expansion_radii(self): + """Return an array of radii associated with the (interleaved) + expansion centers. + + ``coord_t [ncenters]`` + """ + with cl.CommandQueue(self.cl_context) as queue: + from pytential.qbx.utils import get_interleaved_radii + return get_interleaved_radii(queue, self.lpot_source) # }}} @@ -435,22 +449,6 @@ class QBXFMMGeometryData(object): target_discr_starts=target_discr_starts, ntargets=ntargets).with_queue(None) - # FIXME unused - def xxx_target_radii(self): - """Shape: ``[ntargets]`` - - A list of target radii for FMM tree construction. In this list, the QBX - centers have nonzero radii, and all other targets have radius zero. - """ - - tgt_info = self.target_info() - - with cl.CommandQueue(self.cl_context) as queue: - target_radii = cl.array.zeros(queue, tgt_info.ntargets, self.coord_dtype) - # target_radii[:self.ncenters] = center_info.radii - - return target_radii.with_queue(None) - def target_side_preferences(self): """Return one big array combining all the data from the *side* part of :attr:`TargetInfo.target_discrs_and_qbx_sides`. @@ -492,6 +490,12 @@ class QBXFMMGeometryData(object): nsources = lpot_src.fine_density_discr.nnodes nparticles = nsources + target_info.ntargets + target_radii = None + if self.lpot_source.expansion_disks_in_tree_have_extent: + target_radii = cl.array.zeros(queue, target_info.ntargets, + self.coord_dtype) + target_radii[:self.ncenters] = self.expansion_radii() + refine_weights = cl.array.zeros(queue, nparticles, dtype=np.int32) refine_weights[:nsources] = 1 refine_weights.finish() @@ -506,13 +510,16 @@ class QBXFMMGeometryData(object): # to avoid having too many disks overlap more than one box. # # FIXME: Should investigate this further. + tree, _ = code_getter.build_tree(queue, particles=lpot_src.fine_density_discr.nodes(), targets=target_info.targets, + target_radii=target_radii, max_leaf_refine_weight=256, refine_weights=refine_weights, debug=self.debug, - kind="adaptive-level-restricted") + stick_out_factor=0, + kind="adaptive") return tree @@ -530,6 +537,9 @@ class QBXFMMGeometryData(object): trav, _ = self.code_getter.build_traversal(queue, self.tree(), debug=self.debug) + if self.lpot_source.expansion_disks_in_tree_have_extent: + trav = trav.merge_close_lists(queue) + return trav @memoize_method diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index b29f11ab..02a19e00 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -67,6 +67,76 @@ QBX_TREE_MAKO_DEFS = r"""//CL:mako// # }}} +# {{{ interleaver kernel + +@memoize +def get_interleaver_kernel(dtype): + # NOTE: Returned kernel needs dstlen or dst parameter + from pymbolic import var + knl = lp.make_kernel( + "[srclen,dstlen] -> {[i]: 0<=i {[i]: 0<=i`_. + """ + + def __init__(self, domain_n_exprs, ne, + interfaces, use_l2_weighting=None): + """ + :attr interfaces: a tuple of tuples + ``(outer_domain, inner_domain, interface_id)``, + where *outer_domain* and *inner_domain* are indices into + *domain_k_names*, + and *interface_id* is a symbolic name for the discretization of the + interface. 'outer' designates the side of the interface to which + the normal points. + :attr domain_n_exprs: a tuple of variable names of the Helmholtz + parameter *k*, to be used inside each part of the source geometry. + May also be a tuple of strings, which will be transformed into + variable references of the corresponding names. + :attr beta: A symbolic expression for the wave number in the :math:`z` + direction. May be a string, which will be interpreted as a variable + name. + """ + + self.interfaces = interfaces + + ne = sym.var(ne) + self.ne = sym.cse(ne, "ne") + + self.domain_n_exprs = [ + sym.var(n_expr) + for idom, n_expr in enumerate(domain_n_exprs)] + del domain_n_exprs + + import pymbolic.primitives as p + + def upper_half_square_root(x): + return p.If( + p.Comparison( + (x**0.5).a.imag, + "<", 0), + 1j*(-x)**0.5, + x**0.5) + + self.domain_K_exprs = [ + sym.cse( + upper_half_square_root(n_expr**2-ne**2), + "K%d" % i) + for i, n_expr in enumerate(self.domain_n_exprs)] + + from sumpy.kernel import HelmholtzKernel + self.kernel = HelmholtzKernel(2, allow_evanescent=True) + + def make_unknown(self, name): + return sym.make_sym_vector(name, len(self.unknown_index)) + + @property + @memoize_method + def unknown_index(self): + result = {} + for i in range(len(self.interfaces)): + for current in ["j", "m"]: + for cur_dir in ["z", "t"]: + result[(current + cur_dir, i)] = len(result) + + return result + + def tangent(self, where): + return sym.cse( + (sym.pseudoscalar(2, 1, where) + / sym.area_element(2, 1, where)), + "tangent") + + def S(self, dom_idx, density, qbx_forced_limit=+1): # noqa + return sym.S( + self.kernel, density, + k=self.domain_K_exprs[dom_idx], + qbx_forced_limit=qbx_forced_limit) + + def D(self, dom_idx, density, qbx_forced_limit="avg"): # noqa + return sym.D( + self.kernel, density, + k=self.domain_K_exprs[dom_idx], + qbx_forced_limit=qbx_forced_limit) + + def T(self, where, dom_idx, density): # noqa + return sym.int_g_dsource( + 2, + self.tangent(where), + self.kernel, + density, + k=self.domain_K_exprs[dom_idx], + # ???? + qbx_forced_limit=1).xproject(0) + + def representation(self, domain_idx, unknown): + """"Return a tuple of vectors [Ex, Ey, Ez] and [Hx, Hy, Hz] + representing the solution to Maxwell's equation on domain + *domain_idx*. + """ + result = np.zeros(4, dtype=object) + + unk_idx = self.unknown_index + k_v = self.k_vacuum + + for i in range(len(self.interfaces)): + dom0i, dom1i, where = self.interfaces[i] + + if domain_idx == dom0i: + domi = dom0i + elif domain_idx == dom1i: + domi = dom1i + else: + # Interface does not border the requested domain + continue + + beta = self.ne*k_v + k = sym.cse(k_v * self.domain_n_exprs[domi], + "k%d" % domi) + + jt = unknown[unk_idx["jz", i]] + jz = unknown[unk_idx["jt", i]] + mt = unknown[unk_idx["mz", i]] + mz = unknown[unk_idx["mt", i]] + + def diff_axis(iaxis, operand): + v = np.array(2, dtype=np.float64) + v[iaxis] = 1 + d = sym.Derivative() + return d.resolve( + (sym.MultiVector(v).scalar_product(d.dnabla(2))) + * d(operand)) + + from functools import partial + dx = partial(diff_axis, 1) + dy = partial(diff_axis, 1) + + tangent = self.tangent(where) + tau1 = tangent[0] + tau2 = tangent[1] + + S = self.S # noqa + D = self.D # noqa + T = self.T # noqa + + # Ex + result[0] += ( + -1/(1j*k_v) * dx(T(where, domi, jt)) + + beta/k_v * dx(S(domi, jz)) + + k**2/(1j*k_v) * S(domi, jt * tau1) + - dy(S(domi, mz)) + + 1j*beta * S(domi, mt*tau2) + ) + + # Ey + result[1] += ( + - 1/(1j*k_v) * dy(T(where, domi, jt)) + + beta/k_v * dy(S(domi, jz)) + + k**2/(1j*k_v) * S(domi, jt * tau2) + + dx(S(domi, mz)) + - 1j*beta*S(domi, mt * tau1) + ) + + # Ez + result[2] += ( + - beta/k_v * T(where, domi, jt) + + (k**2 - beta**2)/(1j*k_v) * S(domi, jz) + + D(domi, mt) + ) + + # Hx + result[3] += ( + 1/(1j*k_v) * dx(T(where, domi, mt)) + - beta/k_v * dx(S(domi, mz)) + - k**2/(1j*k_v) * S(domi, mt*tau1) + - k**2/k_v**2 * dy(S(domi, jz)) + + 1j * beta * k**2/k_v**2 * S(domi, jt*tau2) + ) + + # Hy + result[4] += ( + 1/(1j*k_v) * dy(T(where, domi, mt)) + - beta/k_v * dy(S(domi, mz)) + - k**2/(1j*k_v) * S(domi, mt * tau2) + + k**2/k_v**2 * dx(S(domi, jz)) + - 1j*beta * k**2/k_v**2 * S(domi, jt*tau1) + ) + + # Hz + result[5] += ( + beta/k_v * T(where, domi, mt) + - (k**2 - beta**2)/(1j*k_v) * S(domi, mz) + + k**2/k_v**2 * D(domi, jt) + ) + + return result + + def operator(self, unknown): + result = np.zeros(4*len(self.interfaces), dtype=object) + + unk_idx = self.unknown_index + + for i in range(len(self.interfaces)): + idx_jt = unk_idx["jz", i] + idx_jz = unk_idx["jt", i] + idx_mt = unk_idx["mz", i] + idx_mz = unk_idx["mt", i] + + phi1 = unknown[idx_jt] + phi2 = unknown[idx_jz] + phi3 = unknown[idx_mt] + phi4 = unknown[idx_mz] + + ne = self.ne + + dom0i, dom1i, where = self.interfaces[i] + + tangent = self.tangent(where) + normal = sym.cse( + sym.normal(2, 1, where), + "normal") + + S = self.S # noqa + D = self.D # noqa + T = self.T # noqa + + def Tt(where, dom, density): # noqa + return sym.tangential_derivative( + 2, T(where, dom, density)).xproject(0) + + def Sn(dom, density): # noqa + return sym.normal_derivative( + 2, + S(dom, density, + qbx_forced_limit="avg")) + + def St(dom, density): # noqa + return sym.tangential_derivative(2, S(dom, density)).xproject(0) + + n0 = self.domain_n_exprs[dom0i] + n1 = self.domain_n_exprs[dom1i] + + a11 = sym.cse(n0**2 * D(dom0i, phi1) - n1**2 * D(dom1i, phi1), "a11") + a22 = sym.cse(-n0**2 * Sn(dom0i, phi2) + n1**2 * Sn(dom1i, phi2), "a22") + a33 = sym.cse(D(dom0i, phi3)-D(dom1i, phi3), "a33") + a44 = sym.cse(-Sn(dom0i, phi4) + Sn(dom1i, phi4), "a44") + + a21 = sym.cse(-1j * ne * ( + n0**2 * tangent.scalar_product( + S(dom0i, normal * phi1)) + - n1**2 * tangent.scalar_product( + S(dom1i, normal * phi1))), "a21") + + a43 = sym.cse(-1j * ne * ( + tangent.scalar_product( + S(dom0i, normal * phi3)) + - tangent.scalar_product( + S(dom1i, normal * phi3))), "a43") + + a13 = +1*sym.cse( + ne*(T(where, dom0i, phi3) - T(where, dom1i, phi3)), "a13") + a31 = -1*sym.cse( + ne*(T(where, dom0i, phi1) - T(where, dom1i, phi1)), "a31") + + a24 = +1*sym.cse(ne*(St(dom0i, phi4) - St(dom1i, phi4)), "a24") + a42 = -1*sym.cse(ne*(St(dom0i, phi2) - St(dom1i, phi2)), "a42") + + a14 = sym.cse(1j*( + (n0**2 - ne**2) * S(dom0i, phi4) + - (n1**2 - ne**2) * S(dom1i, phi4) + ), "a14") + a32 = -sym.cse(1j*( + (n0**2 - ne**2) * S(dom0i, phi2) + - (n1**2 - ne**2) * S(dom1i, phi2) + ), "a32") + + def a23_expr(phi): + return ( + 1j * (Tt(where, dom0i, phi) - Tt(where, dom1i, phi)) + - 1j * ( + n0**2 * tangent.scalar_product( +<<<<<<< HEAD + S(dom0i, normal * phi)) + - n1**2 * tangent.scalar_product( + S(dom1i, normal * phi)))) +======= + S(0, tangent * phi)) + - n1**2 * tangent.scalar_product( + S(1, tangent * phi)))) +>>>>>>> 9af8f2bde0070ee8488c6bc58c8ece5114004a71 + + a23 = +1*sym.cse(a23_expr(phi3), "a23") + a41 = -1*sym.cse(a23_expr(phi1), "a41") + + d1 = (n0**2 + n1**2)/2 * phi1 + d2 = (n0**2 + n1**2)/2 * phi2 + d3 = phi3 + d4 = phi4 + + result[idx_jt] += d1 + a11 + 000 + a13 + a14 + result[idx_jz] += d2 + a21 + a22 + a23 + a24 + result[idx_mt] += d3 + a31 + a32 + a33 + 0 + result[idx_mz] += d4 + a41 + a42 + a43 + a44 + + # TODO: L2 weighting + # TODO: Add representation contributions to other boundaries + # abutting the domain + return result + +# }}} + + +# {{{ old-style waveguide + +class Dielectric2DBoundaryOperatorBase(L2WeightedPDEOperator): + r""" + Solves the following system of BVPs on :math:`\mathbb{R}^2`, in which + a disjoint family of domains :math:`\Omega_i` is embedded: + + .. math:: + + \triangle E + (k_i^2-\beta^2) E = 0\quad \text{on $\Omega_i$}\\ + \triangle H + (k_i^2-\beta^2) H = 0\quad \text{on $\Omega_i$}\\ + [H] = 0 \text{ on $\partial \Omega_i$},\quad + [E] = 0 \text{ on $\partial \Omega_i$}\\ + \left[ \frac{k_0}{k^2-\beta^2} \partial_{\hat n}H\right] = 0 + \quad\text{on $\partial \Omega_i$},\quad\\ + \left[ \frac{k_0}{k^2-\beta^2} \partial_{\hat n}E\right] = 0 + \quad\text{on $\partial \Omega_i$} + + :math:`E` and :math:`H` are assumed to be of the form + + .. math:: + + E(x,y,z,t)=E(x,y)e^{i(\beta z-\omega t) + H(x,y,z,t)=H(x,y)e^{i(\beta z-\omega t) + + where :math:`[\cdot]` denotes the jump across an interface, and :math:`k` + (without an index) denotes the value of :math:`k` on either side of the + interface, for the purpose of computing the jump. :math:`\hat n` denotes + the unit normal of the interface. + + .. automethod:: make_unknown + .. automethod:: representation_outer + .. automethod:: representation_inner + .. automethod:: operator + """ + + field_kind_e = 0 + field_kind_h = 1 + field_kinds = [field_kind_e, field_kind_h] + + side_in = 0 + side_out = 1 + sides = [side_in, side_out] + side_to_sign = { + side_in: -1, + side_out: 1, + } + + dir_none = 0 + dir_normal = 1 + dir_tangential = 2 + + BCTermDescriptor = namedtuple("BCDescriptor", + "i_interface direction field_kind coeff_inner coeff_outer".split()) + + # {{{ constructor + + def __init__(self, mode, k_vacuum, domain_k_exprs, beta, + interfaces, use_l2_weighting=None): + """ + :attr mode: one of 'te', 'tm', 'tem' + :attr k_vacuum: A symbolic expression for the wave number in vacuum. + May be a string, which will be interpreted as a variable name. + :attr interfaces: a tuple of tuples + ``(outer_domain, inner_domain, interface_id)``, + where *outer_domain* and *inner_domain* are indices into + *domain_k_names*, + and *interface_id* is a symbolic name for the discretization of the + interface. 'outer' designates the side of the interface to which + the normal points. + :attr domain_k_exprs: a tuple of variable names of the Helmholtz + parameter *k*, to be used inside each part of the source geometry. + May also be a tuple of strings, which will be transformed into + variable references of the corresponding names. + :attr beta: A symbolic expression for the wave number in the :math:`z` + direction. May be a string, which will be interpreted as a variable + name. + """ + + if use_l2_weighting is None: + use_l2_weighting = False + + super(Dielectric2DBoundaryOperatorBase, self).__init__( + use_l2_weighting=use_l2_weighting) + + if mode == "te": + self.ez_enabled = False + self.hz_enabled = True + elif mode == "tm": + self.ez_enabled = True + self.hz_enabled = False + elif mode == "tem": + self.ez_enabled = True + self.hz_enabled = True + else: + raise ValueError("invalid mode '%s'" % mode) + + self.interfaces = interfaces + + fk_e = self.field_kind_e + fk_h = self.field_kind_h + + dir_none = self.dir_none + dir_normal = self.dir_normal + dir_tangential = self.dir_tangential + + if isinstance(beta, str): + beta = sym.var(beta) + beta = sym.cse(beta, "beta") + + if isinstance(k_vacuum, str): + k_vacuum = sym.var(k_vacuum) + k_vacuum = sym.cse(k_vacuum, "k_vac") + + self.domain_k_exprs = [ + sym.var(k_expr) + if isinstance(k_expr, str) + else sym.cse(k_expr, "k%d" % idom) + for idom, k_expr in enumerate(domain_k_exprs)] + del domain_k_exprs + + # Note the case of k/K! + # "K" is the 2D Helmholtz parameter. + # "k" is the 3D Helmholtz parameter. + + self.domain_K_exprs = [ + sym.cse((k_expr**2-beta**2)**0.5, "K%d" % i) + for i, k_expr in enumerate(self.domain_k_exprs)] + + from sumpy.kernel import HelmholtzKernel + self.kernel = HelmholtzKernel(2, allow_evanescent=True) + + # {{{ build bc list + + # list of tuples, where each tuple consists of BCTermDescriptor instances + + all_bcs = [] + for i_interface, (outer_domain, inner_domain, _) in ( + enumerate(self.interfaces)): + k_outer = self.domain_k_exprs[outer_domain] + k_inner = self.domain_k_exprs[inner_domain] + + all_bcs += [ + ( # [E] = 0 + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_none, + field_kind=fk_e, + coeff_outer=1, + coeff_inner=-1), + ), + ( # [H] = 0 + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_none, + field_kind=fk_h, + coeff_outer=1, + coeff_inner=-1), + ), + ( + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_tangential, + field_kind=fk_e, + coeff_outer=beta/(k_outer**2-beta**2), + coeff_inner=-beta/(k_inner**2-beta**2)), + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_normal, + field_kind=fk_h, + coeff_outer=sym.cse(-k_vacuum/(k_outer**2-beta**2)), + coeff_inner=sym.cse(k_vacuum/(k_inner**2-beta**2))), + ), + ( + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_tangential, + field_kind=fk_h, + coeff_outer=beta/(k_outer**2-beta**2), + coeff_inner=-beta/(k_inner**2-beta**2)), + self.BCTermDescriptor( + i_interface=i_interface, + direction=dir_normal, + field_kind=fk_e, + coeff_outer=sym.cse( + (k_outer**2/k_vacuum)/(k_outer**2-beta**2)), + coeff_inner=sym.cse( + -(k_inner**2/k_vacuum) + / (k_inner**2-beta**2))) + ), + ] + + del k_outer + del k_inner + + self.bcs = [] + for bc in all_bcs: + any_significant_e = any( + term.field_kind == fk_e + and term.direction in [dir_normal, dir_none] + for term in bc) + any_significant_h = any( + term.field_kind == fk_h + and term.direction in [dir_normal, dir_none] + for term in bc) + is_necessary = ( + (self.ez_enabled and any_significant_e) + or + (self.hz_enabled and any_significant_h)) + + # Only keep tangential modes for TEM. Otherwise, + # no jump in H already implies jump condition on + # tangential derivative. + is_tem = self.ez_enabled and self.hz_enabled + terms = tuple( + term + for term in bc + if term.direction != dir_tangential + or is_tem) + + if is_necessary: + self.bcs.append(terms) + + assert (len(all_bcs) + * (int(self.ez_enabled) + int(self.hz_enabled)) // 2 + == len(self.bcs)) + + # }}} + + # }}} + + def is_field_present(self, field_kind): + return ( + (field_kind == self.field_kind_e and self.ez_enabled) + or + (field_kind == self.field_kind_h and self.hz_enabled)) + + def make_unknown(self, name): + num_densities = ( + 2 + * (int(self.ez_enabled) + int(self.hz_enabled)) + * len(self.interfaces)) + + assert num_densities == len(self.bcs) + + return sym.make_sym_vector(name, num_densities) + + def bc_term_to_operator_contrib(self, term, side, raw_potential_op, + density, discrete): + potential_op = raw_potential_op + + side_sign = self.side_to_sign[side] + + domain_outer, domain_inner, interface_id = \ + self.interfaces[term.i_interface] + if side == self.side_in: + K_expr = self.domain_K_exprs[domain_inner] # noqa + bc_coeff = term.coeff_inner + elif side == self.side_out: + K_expr = self.domain_K_exprs[domain_outer] # noqa + bc_coeff = term.coeff_outer + else: + raise ValueError("invalid value of 'side'") + + potential_op = potential_op( + self.kernel, density, source=interface_id, + k=K_expr) + + if term.direction == self.dir_none: + if raw_potential_op is sym.S: + jump_term = 0 + elif raw_potential_op is sym.D: + jump_term = (side_sign*0.5) * discrete + else: + assert False, raw_potential_op + elif term.direction == self.dir_normal: + potential_op = sym.normal_derivative( + potential_op, interface_id) + + if raw_potential_op is sym.S: + # S' + jump_term = (-side_sign*0.5) * discrete + elif raw_potential_op is sym.D: + jump_term = 0 + else: + assert False, raw_potential_op + + elif term.direction == self.dir_tangential: + potential_op = sym.tangential_derivative( + raw_potential_op( + self.kernel, density, source=interface_id, + k=K_expr, qbx_forced_limit=side_sign), + interface_id).a.as_scalar() + + # Some of these may have jumps, but QBX does the dirty + # work here by directly computing the limit. + jump_term = 0 + + else: + raise ValueError("invalid direction") + + potential_op = ( + jump_term + + self.get_sqrt_weight(interface_id)*potential_op) + + del jump_term + + contrib = bc_coeff * potential_op + + if (raw_potential_op is sym.D + and term.direction == self.dir_normal): + # FIXME The hypersingular part should perhaps be + # treated specially to avoid cancellation. + pass + + return contrib + + +# {{{ single-layer representation + +class DielectricSRep2DBoundaryOperator(Dielectric2DBoundaryOperatorBase): + def _structured_unknown(self, unknown, with_l2_weights): + """ + :arg with_l2_weights: If True, return the 'bare' unknowns + that do not have the :math:`L^2` weights divided out. + Note: Those unknowns should *not* be interpreted as + point values of a density. + :returns: an array of unknowns, with the following index axes: + ``[side, field_kind, i_interface]``, where + ``side`` is o for the outside part and i for the interior part, + ``field_kind`` is 0 for the E-field and 1 for the H-field part, + ``i_interface`` is the number of the enclosed domain, starting from 0. + """ + result = np.zeros((2, 2, len(self.interfaces)), dtype=np.object) + + i_unknown = 0 + for side in self.sides: + for field_kind in self.field_kinds: + for i_interface in range(len(self.interfaces)): + + if self.is_field_present(field_kind): + dens = unknown[i_unknown] + i_unknown += 1 + else: + dens = 0 + + _, _, interface_id = self.interfaces[i_interface] + + if not with_l2_weights: + dens = sym.cse( + dens/self.get_sqrt_weight(interface_id), + "dens_{side}_{field}_{dom}".format( + side={ + self.side_out: "o", + self.side_in: "i"} + [side], + field={ + self.field_kind_e: "E", + self.field_kind_h: "H" + } + [field_kind], + dom=i_interface)) + + result[side, field_kind, i_interface] = dens + + assert i_unknown == len(unknown) + return result + + def representation(self, unknown, i_domain): + """ + :return: a symbolic expression for the representation of the PDE solution + in domain number *i_domain*. + """ + unk = self._structured_unknown(unknown, with_l2_weights=False) + + result = [] + + for field_kind in self.field_kinds: + if not self.is_field_present(field_kind): + continue + + field_result = 0 + for i_interface, (i_domain_outer, i_domain_inner, interface_id) in ( + enumerate(self.interfaces)): + if i_domain_outer == i_domain: + side = self.side_out + elif i_domain_inner == i_domain: + side = self.side_in + else: + continue + + my_unk = unk[side, field_kind, i_interface] + if my_unk: + field_result += sym.S( + self.kernel, + my_unk, + source=interface_id, + k=self.domain_K_exprs[i_domain]) + + result.append(field_result) + + from pytools.obj_array import make_obj_array + return make_obj_array(result) + + def operator(self, unknown): + density_unk = self._structured_unknown(unknown, with_l2_weights=False) + discrete_unk = self._structured_unknown(unknown, with_l2_weights=True) + + result = [] + for bc in self.bcs: + op = 0 + + for side in self.sides: + for term in bc: + unk_index = (side, term.field_kind, term.i_interface) + density = density_unk[unk_index] + discrete = discrete_unk[unk_index] + + op += self.bc_term_to_operator_contrib( + term, side, sym.S, density, discrete) + + result.append(op) + + return np.array(result, dtype=np.object) + +# }}} + + +# {{{ single + double layer representation + +class DielectricSDRep2DBoundaryOperator(Dielectric2DBoundaryOperatorBase): + pot_kind_S = 0 + pot_kind_D = 1 + pot_kinds = [pot_kind_S, pot_kind_D] + potential_ops = { + pot_kind_S: sym.S, + pot_kind_D: sym.D, + } + + def __init__(self, mode, k_vacuum, domain_k_exprs, beta, + interfaces, use_l2_weighting=None): + + super(DielectricSDRep2DBoundaryOperator, self).__init__( + mode, k_vacuum, domain_k_exprs, beta, + interfaces, use_l2_weighting=use_l2_weighting) + + side_in = self.side_in + side_out = self.side_out + + def find_normal_derivative_bc_coeff(field_kind, i_interface, side): + result = 0 + for bc in self.bcs: + for term in bc: + if term.field_kind != field_kind: + continue + if term.i_interface != i_interface: + continue + if term.direction != self.dir_normal: + continue + + if side == side_in: + result += term.coeff_inner + elif side == side_out: + result += term.coeff_outer + else: + raise ValueError("invalid side") + + return result + + self.density_coeffs = np.zeros( + (len(self.pot_kinds), len(self.field_kinds), + len(self.interfaces), len(self.sides)), + dtype=np.object) + for field_kind in self.field_kinds: + for i_interface in range(len(self.interfaces)): + self.density_coeffs[ + self.pot_kind_S, field_kind, i_interface, side_in] = 1 + self.density_coeffs[ + self.pot_kind_S, field_kind, i_interface, side_out] = 1 + + # These need to satisfy + # + # [dens_coeff_D * bc_coeff * dn D] + # = dens_coeff_D_out * bc_coeff_out * (dn D) + # + dens_coeff_D_in * bc_coeff_in * dn D + # = 0 + # + # (because dn D is hypersingular, which we'd like to cancel out) + # + # NB: bc_coeff_{in,out} already contain the signs to realize + # the subtraction for the jump. (So the "+" above is as it + # should be.) + + dens_coeff_D_in = find_normal_derivative_bc_coeff( # noqa + field_kind, i_interface, side_out) + dens_coeff_D_out = - find_normal_derivative_bc_coeff( # noqa + field_kind, i_interface, side_in) + + self.density_coeffs[ + self.pot_kind_D, field_kind, i_interface, side_in] \ + = dens_coeff_D_in + self.density_coeffs[ + self.pot_kind_D, field_kind, i_interface, side_out] \ + = dens_coeff_D_out + + def _structured_unknown(self, unknown, with_l2_weights): + """ + :arg with_l2_weights: If True, return the 'bare' unknowns + that do not have the :math:`L^2` weights divided out. + Note: Those unknowns should *not* be interpreted as + point values of a density. + :returns: an array of unknowns, with the following index axes: + ``[pot_kind, field_kind, i_interface]``, where + ``pot_kind`` is 0 for the single-layer part and 1 for the double-layer + part, + ``field_kind`` is 0 for the E-field and 1 for the H-field part, + ``i_interface`` is the number of the enclosed domain, starting from 0. + """ + result = np.zeros((2, 2, len(self.interfaces)), dtype=np.object) + + i_unknown = 0 + for pot_kind in self.pot_kinds: + for field_kind in self.field_kinds: + for i_interface in range(len(self.interfaces)): + + if self.is_field_present(field_kind): + dens = unknown[i_unknown] + i_unknown += 1 + else: + dens = 0 + + _, _, interface_id = self.interfaces[i_interface] + + if not with_l2_weights: + dens = sym.cse( + dens/self.get_sqrt_weight(interface_id), + "dens_{pot}_{field}_{intf}".format( + pot={0: "S", 1: "D"}[pot_kind], + field={ + self.field_kind_e: "E", + self.field_kind_h: "H" + } + [field_kind], + intf=i_interface)) + + result[pot_kind, field_kind, i_interface] = dens + + assert i_unknown == len(unknown) + return result + + def representation(self, unknown, i_domain): + """ + :return: a symbolic expression for the representation of the PDE solution + in domain number *i_domain*. + """ + unk = self._structured_unknown(unknown, with_l2_weights=False) + + result = [] + + for field_kind in self.field_kinds: + if not self.is_field_present(field_kind): + continue + + field_result = 0 + for pot_kind in self.pot_kinds: + for i_interface, (i_domain_outer, i_domain_inner, interface_id) in ( + enumerate(self.interfaces)): + if i_domain_outer == i_domain: + side = self.side_out + elif i_domain_inner == i_domain: + side = self.side_in + else: + continue + + my_unk = unk[pot_kind, field_kind, i_interface] + if my_unk: + field_result += ( + self.density_coeffs[ + pot_kind, field_kind, i_interface, side] + * self.potential_ops[pot_kind]( + self.kernel, + my_unk, + source=interface_id, + k=self.domain_K_exprs[i_domain] + )) + + result.append(field_result) + + from pytools.obj_array import make_obj_array + return make_obj_array(result) + + def operator(self, unknown): + density_unk = self._structured_unknown(unknown, with_l2_weights=False) + discrete_unk = self._structured_unknown(unknown, with_l2_weights=True) + + result = [] + for bc in self.bcs: + op = 0 + + for pot_kind in self.pot_kinds: + for term in bc: + + for side in self.sides: + raw_potential_op = \ + self.potential_ops[pot_kind] + + unk_index = (pot_kind, term.field_kind, term.i_interface) + density = density_unk[unk_index] + discrete = discrete_unk[unk_index] + + op += ( + self.density_coeffs[ + pot_kind, term.field_kind, term.i_interface, + side] + * self.bc_term_to_operator_contrib( + term, side, raw_potential_op, density, discrete) + ) + + result.append(op) + + return np.array(result, dtype=np.object) + +# }}} + +# }}} + +# vim: foldmethod=marker -- GitLab