From 4b6860f7bbc5dc0be031096fb7f2f3b5c0471fd5 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 12 Aug 2017 19:12:19 -0500 Subject: [PATCH 1/4] Rename discretizations flying around the QBX source for (some) consistency --- pytential/qbx/__init__.py | 44 +++++++++++++++++++----------------- pytential/qbx/geometry.py | 8 +++---- pytential/qbx/refinement.py | 18 +++++++-------- pytential/qbx/utils.py | 27 +++++++++++----------- pytential/source.py | 6 ++--- pytential/symbolic/matrix.py | 2 +- pytential/unregularized.py | 13 ++++++----- test/test_global_qbx.py | 5 ++-- 8 files changed, 64 insertions(+), 59 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 50fb8e58..d708a053 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -67,7 +67,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): qbx_order=None, fmm_order=None, fmm_level_to_order=None, - base_resampler=None, + to_refined_connection=None, expansion_factory=None, target_association_tolerance=_not_provided, @@ -83,10 +83,10 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): """ :arg fine_order: The total degree to which the (upsampled) underlying quadrature is exact. - :arg base_resampler: A connection used for resampling from + :arg to_refined_connection: A connection used for resampling from *density_discr* the fine density discretization. It is assumed that the fine density discretization given by - *base_resampler.to_discr* is *not* already upsampled. May + *to_refined_connection.to_discr* is *not* already upsampled. May be *None*. :arg fmm_order: `False` for direct calculation. ``None`` will set a reasonable(-ish?) default. @@ -146,7 +146,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): self.fmm_backend = fmm_backend # Default values are lazily provided if these are None - self._base_resampler = base_resampler + self._to_refined_connection = to_refined_connection if expansion_factory is None: from sumpy.expansion import DefaultExpansionFactory @@ -166,7 +166,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): fine_order=None, qbx_order=None, fmm_level_to_order=None, - base_resampler=None, + to_refined_connection=None, target_association_tolerance=_not_provided, _expansions_in_tree_have_extent=_not_provided, _expansion_stick_out_factor=_not_provided, @@ -210,7 +210,8 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): fmm_level_to_order=( fmm_level_to_order or self.fmm_level_to_order), target_association_tolerance=target_association_tolerance, - base_resampler=base_resampler or self._base_resampler, + to_refined_connection=( + to_refined_connection or self._to_refined_connection), debug=( # False is a valid value here @@ -237,25 +238,23 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): # }}} @property - def base_fine_density_discr(self): + def refined_interp_density_discr(self): """The refined, interpolation-focused density discretization (no oversampling). """ - # FIXME: Maybe rename refined_interp_density_discr - return (self._base_resampler.to_discr - if self._base_resampler is not None + return (self._to_refined_connection.to_discr + if self._to_refined_connection is not None else self.density_discr) @property @memoize_method - def fine_density_discr(self): + def refined_ovsmp_quad_density_discr(self): """The refined, quadrature-focused density discretization (with upsampling). """ - # FIXME: Maybe rename refined_quad_density_discr from meshmode.discretization.poly_element import ( QuadratureSimplexGroupFactory) return Discretization( - self.density_discr.cl_context, self.base_fine_density_discr.mesh, + self.density_discr.cl_context, self.refined_interp_density_discr.mesh, QuadratureSimplexGroupFactory(self.fine_order), self.real_dtype) @@ -266,10 +265,12 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): make_same_mesh_connection, ChainedDiscretizationConnection) conn = make_same_mesh_connection( - self.fine_density_discr, self.base_fine_density_discr) + self.refined_ovsmp_quad_density_discr, + self.refined_interp_density_discr) - if self._base_resampler is not None: - return ChainedDiscretizationConnection([self._base_resampler, conn]) + if self._to_refined_connection is not None: + return ChainedDiscretizationConnection( + [self._to_refined_connection, conn]) return conn @@ -329,7 +330,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): @memoize_method def _fine_panel_centers_of_mass(self): import pytential.qbx.utils as utils - return utils.element_centers_of_mass(self.base_fine_density_discr) + return utils.element_centers_of_mass(self.refined_interp_density_discr) @memoize_method def _expansion_radii(self, last_dim_length): @@ -364,7 +365,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): raise NotImplementedError() import pytential.qbx.utils as utils - return utils.panel_sizes(self.base_fine_density_discr, last_dim_length) + return utils.panel_sizes(self.refined_interp_density_discr, last_dim_length) @memoize_method def qbx_fmm_geometry_data(self, target_discrs_and_qbx_sides): @@ -646,7 +647,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): evt, output_for_each_kernel = lpot_applier( queue, target_discr.nodes(), - self.fine_density_discr.nodes(), + self.refined_ovsmp_quad_density_discr.nodes(), utils.get_centers_on_side(self, o.qbx_forced_limit), [strengths], **kernel_args) result.append((o.name, output_for_each_kernel[o.kernel_index])) @@ -659,7 +660,8 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): insn.kernels) evt, output_for_each_kernel = p2p(queue, - target_discr.nodes(), self.fine_density_discr.nodes(), + target_discr.nodes(), + self.refined_ovsmp_quad_density_discr.nodes(), [strengths], **kernel_args) qbx_forced_limit = o.qbx_forced_limit @@ -708,7 +710,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): lpot_applier_on_tgt_subset( queue, targets=target_discr.nodes(), - sources=self.fine_density_discr.nodes(), + sources=self.refined_ovsmp_quad_density_discr.nodes(), centers=geo_data.centers(), strengths=[strengths], qbx_tgt_numbers=qbx_tgt_numbers, diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 96bacdf1..deaf6cf3 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -368,7 +368,7 @@ class QBXFMMGeometryData(object): @property def coord_dtype(self): - return self.lpot_source.fine_density_discr.nodes().dtype + return self.lpot_source.refined_ovsmp_quad_density_discr.nodes().dtype # {{{ centers/radii @@ -480,7 +480,7 @@ class QBXFMMGeometryData(object): target_info = self.target_info() with cl.CommandQueue(self.cl_context) as queue: - nsources = lpot_src.fine_density_discr.nnodes + nsources = lpot_src.refined_ovsmp_quad_density_discr.nnodes nparticles = nsources + target_info.ntargets target_radii = None @@ -505,7 +505,7 @@ class QBXFMMGeometryData(object): # FIXME: Should investigate this further. tree, _ = code_getter.build_tree(queue, - particles=lpot_src.fine_density_discr.nodes(), + particles=lpot_src.refined_ovsmp_quad_density_discr.nodes(), targets=target_info.targets, target_radii=target_radii, max_leaf_refine_weight=32, @@ -806,7 +806,7 @@ class QBXFMMGeometryData(object): with cl.CommandQueue(self.cl_context) as queue: import matplotlib.pyplot as pt from meshmode.discretization.visualization import draw_curve - draw_curve(self.lpot_source.fine_density_discr) + draw_curve(self.lpot_source.refined_ovsmp_quad_density_discr) global_flags = self.global_qbx_flags().get(queue=queue) diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 6c47a104..04eb6487 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -458,7 +458,7 @@ class RefinerNotConvergedWarning(UserWarning): pass -def make_empty_refine_flags(queue, lpot_source, use_base_fine_discr=False): +def make_empty_refine_flags(queue, lpot_source, use_refined_interp_discr=False): """Return an array on the device suitable for use as element refine flags. :arg queue: An instance of :class:`pyopencl.CommandQueue`. @@ -467,8 +467,8 @@ def make_empty_refine_flags(queue, lpot_source, use_base_fine_discr=False): :returns: A :class:`pyopencl.array.Array` suitable for use as refine flags, initialized to zero. """ - discr = (lpot_source.base_fine_density_discr - if use_base_fine_discr + discr = (lpot_source.refined_interp_density_discr + if use_refined_interp_discr else lpot_source.density_discr) result = cl.array.zeros(queue, discr.mesh.nelements, np.int32) result.finish() @@ -578,7 +578,7 @@ def refine_for_global_qbx(lpot_source, wrangler, niter = 0 fine_connections = [] - base_fine_density_discr = lpot_source.density_discr + refined_interp_density_discr = lpot_source.density_discr while must_refine: must_refine = False @@ -594,22 +594,22 @@ def refine_for_global_qbx(lpot_source, wrangler, # Build tree and auxiliary data. # FIXME: The tree should not have to be rebuilt at each iteration. - tree = wrangler.build_tree(lpot_source, use_base_fine_discr=True) + tree = wrangler.build_tree(lpot_source, use_refined_interp_discr=True) peer_lists = wrangler.find_peer_lists(tree) refine_flags = make_empty_refine_flags( - wrangler.queue, lpot_source, use_base_fine_discr=True) + wrangler.queue, lpot_source, use_refined_interp_discr=True) must_refine |= wrangler.check_sufficient_source_quadrature_resolution( lpot_source, tree, peer_lists, refine_flags, debug) if must_refine: conn = wrangler.refine( - base_fine_density_discr, + refined_interp_density_discr, refiner, refine_flags, group_factory, debug) - base_fine_density_discr = conn.to_discr + refined_interp_density_discr = conn.to_discr fine_connections.append(conn) lpot_source = lpot_source.copy( - base_resampler=ChainedDiscretizationConnection( + to_refined_connection=ChainedDiscretizationConnection( fine_connections)) del tree diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index edf92dfe..53b45d52 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -142,12 +142,12 @@ def get_interleaved_radii(queue, lpot_source): class TreeWranglerBase(object): def build_tree(self, lpot_source, targets_list=(), - use_base_fine_discr=False): + use_refined_interp_discr=False): tb = self.code_container.tree_builder() from pytential.qbx.utils import build_tree_with_qbx_metadata return build_tree_with_qbx_metadata( self.queue, tb, lpot_source, targets_list=targets_list, - use_base_fine_discr=use_base_fine_discr) + use_refined_interp_discr=use_refined_interp_discr) def find_peer_lists(self, tree): plf = self.code_container.peer_list_finder() @@ -413,12 +413,13 @@ MAX_REFINE_WEIGHT = 64 def build_tree_with_qbx_metadata( queue, tree_builder, lpot_source, targets_list=(), - use_base_fine_discr=False): + use_refined_interp_discr=False): """Return a :class:`TreeWithQBXMetadata` built from the given layer potential source. This contains particles of four different types: * source particles and panel centers of mass either from - ``lpot_source.density_discr`` or ``lpot_source.base_fine_density_discr`` + ``lpot_source.density_discr`` or + ``lpot_source.refined_interp_density_discr`` * centers from ``lpot_source.centers()`` * targets from ``targets_list``. @@ -429,8 +430,8 @@ def build_tree_with_qbx_metadata( :arg targets_list: A list of :class:`pytential.target.TargetBase` - :arg use_base_fine_discr: If *True*, builds a tree with sources/centers of - mass from ``lpot_source.base_fine_density_discr``. If *False* (default), + :arg use_refined_interp_discr: If *True*, builds a tree with sources/centers of + mass from ``lpot_source.refined_interp_density_discr``. If *False* (default), they are from ``lpot_source.density_discr``. """ # The ordering of particles is as follows: @@ -443,14 +444,14 @@ def build_tree_with_qbx_metadata( sources = ( lpot_source.density_discr.nodes() - if not use_base_fine_discr - else lpot_source.fine_density_discr.nodes()) + if not use_refined_interp_discr + else lpot_source.refined_ovsmp_quad_density_discr.nodes()) centers = get_interleaved_centers(queue, lpot_source) centers_of_mass = ( lpot_source._panel_centers_of_mass() - if not use_base_fine_discr + if not use_refined_interp_discr else lpot_source._fine_panel_centers_of_mass()) targets = (tgt.nodes() for tgt in targets_list) @@ -465,7 +466,7 @@ def build_tree_with_qbx_metadata( nsources = len(sources[0]) ncenters = len(centers[0]) # Each source gets an interior / exterior center. - assert 2 * nsources == ncenters or use_base_fine_discr + assert 2 * nsources == ncenters or use_refined_interp_discr ntargets = sum(tgt.nnodes for tgt in targets_list) # Slices @@ -519,8 +520,8 @@ def build_tree_with_qbx_metadata( del box_to_class # Compute panel => source relation - if use_base_fine_discr: - density_discr = lpot_source.fine_density_discr + if use_refined_interp_discr: + density_discr = lpot_source.refined_ovsmp_quad_density_discr else: density_discr = lpot_source.density_discr @@ -539,7 +540,7 @@ def build_tree_with_qbx_metadata( # Compute panel => center relation qbx_panel_to_center_starts = ( 2 * qbx_panel_to_source_starts - if not use_base_fine_discr + if not use_refined_interp_discr else None) # Transfer all tree attributes. diff --git a/pytential/source.py b/pytential/source.py index 75d0db70..29fafaf1 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -151,7 +151,7 @@ class LayerPotentialSourceBase(PotentialSource): .. rubric:: Discretizations .. attribute:: density_discr - .. attribute:: fine_density_discr + .. attribute:: refined_ovsmp_quad_density_discr .. attribute:: resampler .. method:: with_refinement @@ -267,7 +267,7 @@ class LayerPotentialSourceBase(PotentialSource): import pytential.symbolic.primitives as p from pytential.symbolic.execution import bind with cl.CommandQueue(self.cl_context) as queue: - # fine_density_discr is not guaranteed to be usable for + # refined_ovsmp_quad_density_discr is not guaranteed to be usable for # interpolation/differentiation. Use density_discr to find # area element instead, then upsample that. @@ -277,7 +277,7 @@ class LayerPotentialSourceBase(PotentialSource): p.area_element(self.ambient_dim, self.dim) )(queue)) - qweight = bind(self.fine_density_discr, p.QWeight())(queue) + qweight = bind(self.refined_ovsmp_quad_density_discr, p.QWeight())(queue) return (area_element.with_queue(queue)*qweight).with_queue(None) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 2012b992..727a8d46 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -188,7 +188,7 @@ class MatrixBuilder(EvaluationMapperBase): assert abs(expr.qbx_forced_limit) > 0 _, (mat,) = mat_gen(self.queue, target_discr.nodes(), - source.fine_density_discr.nodes(), + source.refined_ovsmp_quad_density_discr.nodes(), get_centers_on_side(source, expr.qbx_forced_limit), **kernel_args) diff --git a/pytential/unregularized.py b/pytential/unregularized.py index fcd5298e..6616b190 100644 --- a/pytential/unregularized.py +++ b/pytential/unregularized.py @@ -89,7 +89,7 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): self.debug = debug @property - def fine_density_discr(self): + def refined_ovsmp_quad_density_discr(self): return self.density_discr def resampler(self, queue, f): @@ -160,7 +160,8 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): p2p = self.get_p2p(insn.kernels) evt, output_for_each_kernel = p2p(queue, - target_discr.nodes(), self.fine_density_discr.nodes(), + target_discr.nodes(), + self.refined_ovsmp_quad_density_discr.nodes(), [strengths], **kernel_args) result.append((o.name, output_for_each_kernel[o.kernel_index])) @@ -342,11 +343,11 @@ class _FMMGeometryData(object): @property def coord_dtype(self): - return self.lpot_source.fine_density_discr.nodes().dtype + return self.lpot_source.refined_ovsmp_quad_density_discr.nodes().dtype @property def ambient_dim(self): - return self.lpot_source.fine_density_discr.ambient_dim + return self.lpot_source.refined_ovsmp_quad_density_discr.ambient_dim @memoize_method def traversal(self): @@ -369,7 +370,7 @@ class _FMMGeometryData(object): target_info = self.target_info() with cl.CommandQueue(self.cl_context) as queue: - nsources = lpot_src.fine_density_discr.nnodes + nsources = lpot_src.refined_ovsmp_quad_density_discr.nnodes nparticles = nsources + target_info.ntargets refine_weights = cl.array.zeros(queue, nparticles, dtype=np.int32) @@ -379,7 +380,7 @@ class _FMMGeometryData(object): MAX_LEAF_REFINE_WEIGHT = 32 # noqa tree, _ = code_getter.build_tree(queue, - particles=lpot_src.fine_density_discr.nodes(), + particles=lpot_src.refined_ovsmp_quad_density_discr.nodes(), targets=target_info.targets, max_leaf_refine_weight=MAX_LEAF_REFINE_WEIGHT, refine_weights=refine_weights, diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index 12c816f4..fcb172d2 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -108,7 +108,8 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): 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) + fine_discr_nodes = \ + lpot_source.refined_ovsmp_quad_density_discr.nodes().get(queue) int_centers = get_centers_on_side(lpot_source, -1) int_centers = np.array([axis.get(queue) for axis in int_centers]) ext_centers = get_centers_on_side(lpot_source, +1) @@ -176,7 +177,7 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): for i, panel_1 in enumerate(iter_elements(lpot_source.density_discr)): for panel_2 in iter_elements(lpot_source.density_discr): check_disk_undisturbed_by_sources(panel_1, panel_2) - for panel_2 in iter_elements(lpot_source.fine_density_discr): + for panel_2 in iter_elements(lpot_source.refined_ovsmp_quad_density_discr): check_sufficient_quadrature_resolution(panel_1, panel_2) if helmholtz_k is not None: check_panel_size_to_helmholtz_k_ratio(panel_1) -- GitLab From 03ff456af56025e231413b246e2c4270827887da Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 12 Aug 2017 19:29:40 -0500 Subject: [PATCH 2/4] Use refined (not unrefined) density discretization to compute area elements (Closes #66 on Gitlab) --- pytential/qbx/__init__.py | 17 ++++++++++++----- pytential/source.py | 6 ++++-- 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index d708a053..61f500ad 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -245,6 +245,15 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): if self._to_refined_connection is not None else self.density_discr) + @property + @memoize_method + def refined_interp_to_ovsmp_quad_connection(self): + from meshmode.discretization.connection import make_same_mesh_connection + + return make_same_mesh_connection( + self.refined_ovsmp_quad_density_discr, + self.refined_interp_density_discr) + @property @memoize_method def refined_ovsmp_quad_density_discr(self): @@ -261,12 +270,10 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): @property @memoize_method def resampler(self): - from meshmode.discretization.connection import ( - make_same_mesh_connection, ChainedDiscretizationConnection) + from meshmode.discretization.connection import \ + ChainedDiscretizationConnection - conn = make_same_mesh_connection( - self.refined_ovsmp_quad_density_discr, - self.refined_interp_density_discr) + conn = self.refined_interp_to_ovsmp_quad_connection if self._to_refined_connection is not None: return ChainedDiscretizationConnection( diff --git a/pytential/source.py b/pytential/source.py index 29fafaf1..6ec92e32 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -151,6 +151,7 @@ class LayerPotentialSourceBase(PotentialSource): .. rubric:: Discretizations .. attribute:: density_discr + .. attribute:: refined_interp_density_discr .. attribute:: refined_ovsmp_quad_density_discr .. attribute:: resampler .. method:: with_refinement @@ -271,9 +272,10 @@ class LayerPotentialSourceBase(PotentialSource): # interpolation/differentiation. Use density_discr to find # area element instead, then upsample that. - area_element = self.resampler(queue, + area_element = self.refined_interp_to_ovsmp_quad_connection( + queue, bind( - self.density_discr, + self.refined_interp_density_discr, p.area_element(self.ambient_dim, self.dim) )(queue)) -- GitLab From f6f5c98c8282bb8164b6b2fdc43efb5664c84b72 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 13 Aug 2017 02:35:46 -0500 Subject: [PATCH 3/4] Sink weights_and_area_elements() implementation into subclasses, use the old behavior for the unregularized case. --- pytential/qbx/__init__.py | 24 ++++++++++++++++++++++++ pytential/source.py | 23 ----------------------- pytential/unregularized.py | 20 ++++++++++++++++++++ 3 files changed, 44 insertions(+), 23 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 61f500ad..0afb9a22 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -267,6 +267,30 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): QuadratureSimplexGroupFactory(self.fine_order), self.real_dtype) + # {{{ weights and area elements + + @memoize_method + def weights_and_area_elements(self): + import pytential.symbolic.primitives as p + from pytential.symbolic.execution import bind + with cl.CommandQueue(self.cl_context) as queue: + # refined_ovsmp_quad_density_discr is not guaranteed to be usable for + # interpolation/differentiation. Use density_discr to find + # area element instead, then upsample that. + + area_element = self.refined_interp_to_ovsmp_quad_connection( + queue, + bind( + self.refined_interp_density_discr, + p.area_element(self.ambient_dim, self.dim) + )(queue)) + + qweight = bind(self.refined_ovsmp_quad_density_discr, p.QWeight())(queue) + + return (area_element.with_queue(queue)*qweight).with_queue(None) + + # }}} + @property @memoize_method def resampler(self): diff --git a/pytential/source.py b/pytential/source.py index 6ec92e32..91f43db0 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -261,28 +261,5 @@ class LayerPotentialSourceBase(PotentialSource): # }}} - # {{{ weights and area elements - - @memoize_method - def weights_and_area_elements(self): - import pytential.symbolic.primitives as p - from pytential.symbolic.execution import bind - with cl.CommandQueue(self.cl_context) as queue: - # refined_ovsmp_quad_density_discr is not guaranteed to be usable for - # interpolation/differentiation. Use density_discr to find - # area element instead, then upsample that. - - area_element = self.refined_interp_to_ovsmp_quad_connection( - queue, - bind( - self.refined_interp_density_discr, - p.area_element(self.ambient_dim, self.dim) - )(queue)) - - qweight = bind(self.refined_ovsmp_quad_density_discr, p.QWeight())(queue) - - return (area_element.with_queue(queue)*qweight).with_queue(None) - - # }}} # }}} diff --git a/pytential/unregularized.py b/pytential/unregularized.py index 6616b190..16689dfb 100644 --- a/pytential/unregularized.py +++ b/pytential/unregularized.py @@ -88,6 +88,26 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): self.debug = debug + @memoize_method + def weights_and_area_elements(self): + import pytential.symbolic.primitives as p + from pytential.symbolic.execution import bind + with cl.CommandQueue(self.cl_context) as queue: + # refined_ovsmp_quad_density_discr is not guaranteed to be usable for + # interpolation/differentiation. Use density_discr to find + # area element instead, then upsample that. + + area_element = self.resampler( + queue, + bind( + self.density_discr, + p.area_element(self.ambient_dim, self.dim) + )(queue)) + + qweight = bind(self.refined_ovsmp_quad_density_discr, p.QWeight())(queue) + + return (area_element.with_queue(queue)*qweight).with_queue(None) + @property def refined_ovsmp_quad_density_discr(self): return self.density_discr -- GitLab From 79a16a86ad455c8631182f224ff1444515a8c5ab Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 13 Aug 2017 18:59:20 -0500 Subject: [PATCH 4/4] Second round of stage 2 density discr renames --- pytential/qbx/__init__.py | 26 +++++++++++++------------- pytential/qbx/geometry.py | 8 ++++---- pytential/qbx/refinement.py | 16 ++++++++-------- pytential/qbx/utils.py | 26 +++++++++++++------------- pytential/source.py | 4 ++-- pytential/symbolic/matrix.py | 2 +- pytential/unregularized.py | 16 ++++++++-------- test/test_global_qbx.py | 4 ++-- 8 files changed, 51 insertions(+), 51 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 0afb9a22..f0a0a6fb 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -238,7 +238,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): # }}} @property - def refined_interp_density_discr(self): + def stage2_density_discr(self): """The refined, interpolation-focused density discretization (no oversampling). """ return (self._to_refined_connection.to_discr @@ -251,19 +251,19 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): from meshmode.discretization.connection import make_same_mesh_connection return make_same_mesh_connection( - self.refined_ovsmp_quad_density_discr, - self.refined_interp_density_discr) + self.quad_stage2_density_discr, + self.stage2_density_discr) @property @memoize_method - def refined_ovsmp_quad_density_discr(self): + def quad_stage2_density_discr(self): """The refined, quadrature-focused density discretization (with upsampling). """ from meshmode.discretization.poly_element import ( QuadratureSimplexGroupFactory) return Discretization( - self.density_discr.cl_context, self.refined_interp_density_discr.mesh, + self.density_discr.cl_context, self.stage2_density_discr.mesh, QuadratureSimplexGroupFactory(self.fine_order), self.real_dtype) @@ -274,18 +274,18 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): import pytential.symbolic.primitives as p from pytential.symbolic.execution import bind with cl.CommandQueue(self.cl_context) as queue: - # refined_ovsmp_quad_density_discr is not guaranteed to be usable for + # quad_stage2_density_discr is not guaranteed to be usable for # interpolation/differentiation. Use density_discr to find # area element instead, then upsample that. area_element = self.refined_interp_to_ovsmp_quad_connection( queue, bind( - self.refined_interp_density_discr, + self.stage2_density_discr, p.area_element(self.ambient_dim, self.dim) )(queue)) - qweight = bind(self.refined_ovsmp_quad_density_discr, p.QWeight())(queue) + qweight = bind(self.quad_stage2_density_discr, p.QWeight())(queue) return (area_element.with_queue(queue)*qweight).with_queue(None) @@ -361,7 +361,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): @memoize_method def _fine_panel_centers_of_mass(self): import pytential.qbx.utils as utils - return utils.element_centers_of_mass(self.refined_interp_density_discr) + return utils.element_centers_of_mass(self.stage2_density_discr) @memoize_method def _expansion_radii(self, last_dim_length): @@ -396,7 +396,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): raise NotImplementedError() import pytential.qbx.utils as utils - return utils.panel_sizes(self.refined_interp_density_discr, last_dim_length) + return utils.panel_sizes(self.stage2_density_discr, last_dim_length) @memoize_method def qbx_fmm_geometry_data(self, target_discrs_and_qbx_sides): @@ -678,7 +678,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): evt, output_for_each_kernel = lpot_applier( queue, target_discr.nodes(), - self.refined_ovsmp_quad_density_discr.nodes(), + self.quad_stage2_density_discr.nodes(), utils.get_centers_on_side(self, o.qbx_forced_limit), [strengths], **kernel_args) result.append((o.name, output_for_each_kernel[o.kernel_index])) @@ -692,7 +692,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): evt, output_for_each_kernel = p2p(queue, target_discr.nodes(), - self.refined_ovsmp_quad_density_discr.nodes(), + self.quad_stage2_density_discr.nodes(), [strengths], **kernel_args) qbx_forced_limit = o.qbx_forced_limit @@ -741,7 +741,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): lpot_applier_on_tgt_subset( queue, targets=target_discr.nodes(), - sources=self.refined_ovsmp_quad_density_discr.nodes(), + sources=self.quad_stage2_density_discr.nodes(), centers=geo_data.centers(), strengths=[strengths], qbx_tgt_numbers=qbx_tgt_numbers, diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index deaf6cf3..9a2914ae 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -368,7 +368,7 @@ class QBXFMMGeometryData(object): @property def coord_dtype(self): - return self.lpot_source.refined_ovsmp_quad_density_discr.nodes().dtype + return self.lpot_source.quad_stage2_density_discr.nodes().dtype # {{{ centers/radii @@ -480,7 +480,7 @@ class QBXFMMGeometryData(object): target_info = self.target_info() with cl.CommandQueue(self.cl_context) as queue: - nsources = lpot_src.refined_ovsmp_quad_density_discr.nnodes + nsources = lpot_src.quad_stage2_density_discr.nnodes nparticles = nsources + target_info.ntargets target_radii = None @@ -505,7 +505,7 @@ class QBXFMMGeometryData(object): # FIXME: Should investigate this further. tree, _ = code_getter.build_tree(queue, - particles=lpot_src.refined_ovsmp_quad_density_discr.nodes(), + particles=lpot_src.quad_stage2_density_discr.nodes(), targets=target_info.targets, target_radii=target_radii, max_leaf_refine_weight=32, @@ -806,7 +806,7 @@ class QBXFMMGeometryData(object): with cl.CommandQueue(self.cl_context) as queue: import matplotlib.pyplot as pt from meshmode.discretization.visualization import draw_curve - draw_curve(self.lpot_source.refined_ovsmp_quad_density_discr) + draw_curve(self.lpot_source.quad_stage2_density_discr) global_flags = self.global_qbx_flags().get(queue=queue) diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 04eb6487..e1b74563 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -458,7 +458,7 @@ class RefinerNotConvergedWarning(UserWarning): pass -def make_empty_refine_flags(queue, lpot_source, use_refined_interp_discr=False): +def make_empty_refine_flags(queue, lpot_source, use_stage2_discr=False): """Return an array on the device suitable for use as element refine flags. :arg queue: An instance of :class:`pyopencl.CommandQueue`. @@ -467,8 +467,8 @@ def make_empty_refine_flags(queue, lpot_source, use_refined_interp_discr=False): :returns: A :class:`pyopencl.array.Array` suitable for use as refine flags, initialized to zero. """ - discr = (lpot_source.refined_interp_density_discr - if use_refined_interp_discr + discr = (lpot_source.stage2_density_discr + if use_stage2_discr else lpot_source.density_discr) result = cl.array.zeros(queue, discr.mesh.nelements, np.int32) result.finish() @@ -578,7 +578,7 @@ def refine_for_global_qbx(lpot_source, wrangler, niter = 0 fine_connections = [] - refined_interp_density_discr = lpot_source.density_discr + stage2_density_discr = lpot_source.density_discr while must_refine: must_refine = False @@ -594,19 +594,19 @@ def refine_for_global_qbx(lpot_source, wrangler, # Build tree and auxiliary data. # FIXME: The tree should not have to be rebuilt at each iteration. - tree = wrangler.build_tree(lpot_source, use_refined_interp_discr=True) + tree = wrangler.build_tree(lpot_source, use_stage2_discr=True) peer_lists = wrangler.find_peer_lists(tree) refine_flags = make_empty_refine_flags( - wrangler.queue, lpot_source, use_refined_interp_discr=True) + wrangler.queue, lpot_source, use_stage2_discr=True) must_refine |= wrangler.check_sufficient_source_quadrature_resolution( lpot_source, tree, peer_lists, refine_flags, debug) if must_refine: conn = wrangler.refine( - refined_interp_density_discr, + stage2_density_discr, refiner, refine_flags, group_factory, debug) - refined_interp_density_discr = conn.to_discr + stage2_density_discr = conn.to_discr fine_connections.append(conn) lpot_source = lpot_source.copy( to_refined_connection=ChainedDiscretizationConnection( diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index 53b45d52..c27b9a87 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -142,12 +142,12 @@ def get_interleaved_radii(queue, lpot_source): class TreeWranglerBase(object): def build_tree(self, lpot_source, targets_list=(), - use_refined_interp_discr=False): + use_stage2_discr=False): tb = self.code_container.tree_builder() from pytential.qbx.utils import build_tree_with_qbx_metadata return build_tree_with_qbx_metadata( self.queue, tb, lpot_source, targets_list=targets_list, - use_refined_interp_discr=use_refined_interp_discr) + use_stage2_discr=use_stage2_discr) def find_peer_lists(self, tree): plf = self.code_container.peer_list_finder() @@ -413,13 +413,13 @@ MAX_REFINE_WEIGHT = 64 def build_tree_with_qbx_metadata( queue, tree_builder, lpot_source, targets_list=(), - use_refined_interp_discr=False): + use_stage2_discr=False): """Return a :class:`TreeWithQBXMetadata` built from the given layer potential source. This contains particles of four different types: * source particles and panel centers of mass either from ``lpot_source.density_discr`` or - ``lpot_source.refined_interp_density_discr`` + ``lpot_source.stage2_density_discr`` * centers from ``lpot_source.centers()`` * targets from ``targets_list``. @@ -430,8 +430,8 @@ def build_tree_with_qbx_metadata( :arg targets_list: A list of :class:`pytential.target.TargetBase` - :arg use_refined_interp_discr: If *True*, builds a tree with sources/centers of - mass from ``lpot_source.refined_interp_density_discr``. If *False* (default), + :arg use_stage2_discr: If *True*, builds a tree with sources/centers of + mass from ``lpot_source.stage2_density_discr``. If *False* (default), they are from ``lpot_source.density_discr``. """ # The ordering of particles is as follows: @@ -444,14 +444,14 @@ def build_tree_with_qbx_metadata( sources = ( lpot_source.density_discr.nodes() - if not use_refined_interp_discr - else lpot_source.refined_ovsmp_quad_density_discr.nodes()) + if not use_stage2_discr + else lpot_source.quad_stage2_density_discr.nodes()) centers = get_interleaved_centers(queue, lpot_source) centers_of_mass = ( lpot_source._panel_centers_of_mass() - if not use_refined_interp_discr + if not use_stage2_discr else lpot_source._fine_panel_centers_of_mass()) targets = (tgt.nodes() for tgt in targets_list) @@ -466,7 +466,7 @@ def build_tree_with_qbx_metadata( nsources = len(sources[0]) ncenters = len(centers[0]) # Each source gets an interior / exterior center. - assert 2 * nsources == ncenters or use_refined_interp_discr + assert 2 * nsources == ncenters or use_stage2_discr ntargets = sum(tgt.nnodes for tgt in targets_list) # Slices @@ -520,8 +520,8 @@ def build_tree_with_qbx_metadata( del box_to_class # Compute panel => source relation - if use_refined_interp_discr: - density_discr = lpot_source.refined_ovsmp_quad_density_discr + if use_stage2_discr: + density_discr = lpot_source.quad_stage2_density_discr else: density_discr = lpot_source.density_discr @@ -540,7 +540,7 @@ def build_tree_with_qbx_metadata( # Compute panel => center relation qbx_panel_to_center_starts = ( 2 * qbx_panel_to_source_starts - if not use_refined_interp_discr + if not use_stage2_discr else None) # Transfer all tree attributes. diff --git a/pytential/source.py b/pytential/source.py index 91f43db0..f6cabf20 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -151,8 +151,8 @@ class LayerPotentialSourceBase(PotentialSource): .. rubric:: Discretizations .. attribute:: density_discr - .. attribute:: refined_interp_density_discr - .. attribute:: refined_ovsmp_quad_density_discr + .. attribute:: stage2_density_discr + .. attribute:: quad_stage2_density_discr .. attribute:: resampler .. method:: with_refinement diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 727a8d46..9f2f101a 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -188,7 +188,7 @@ class MatrixBuilder(EvaluationMapperBase): assert abs(expr.qbx_forced_limit) > 0 _, (mat,) = mat_gen(self.queue, target_discr.nodes(), - source.refined_ovsmp_quad_density_discr.nodes(), + source.quad_stage2_density_discr.nodes(), get_centers_on_side(source, expr.qbx_forced_limit), **kernel_args) diff --git a/pytential/unregularized.py b/pytential/unregularized.py index 16689dfb..43ad6e59 100644 --- a/pytential/unregularized.py +++ b/pytential/unregularized.py @@ -93,7 +93,7 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): import pytential.symbolic.primitives as p from pytential.symbolic.execution import bind with cl.CommandQueue(self.cl_context) as queue: - # refined_ovsmp_quad_density_discr is not guaranteed to be usable for + # quad_stage2_density_discr is not guaranteed to be usable for # interpolation/differentiation. Use density_discr to find # area element instead, then upsample that. @@ -104,12 +104,12 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): p.area_element(self.ambient_dim, self.dim) )(queue)) - qweight = bind(self.refined_ovsmp_quad_density_discr, p.QWeight())(queue) + qweight = bind(self.quad_stage2_density_discr, p.QWeight())(queue) return (area_element.with_queue(queue)*qweight).with_queue(None) @property - def refined_ovsmp_quad_density_discr(self): + def quad_stage2_density_discr(self): return self.density_discr def resampler(self, queue, f): @@ -181,7 +181,7 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase): evt, output_for_each_kernel = p2p(queue, target_discr.nodes(), - self.refined_ovsmp_quad_density_discr.nodes(), + self.quad_stage2_density_discr.nodes(), [strengths], **kernel_args) result.append((o.name, output_for_each_kernel[o.kernel_index])) @@ -363,11 +363,11 @@ class _FMMGeometryData(object): @property def coord_dtype(self): - return self.lpot_source.refined_ovsmp_quad_density_discr.nodes().dtype + return self.lpot_source.quad_stage2_density_discr.nodes().dtype @property def ambient_dim(self): - return self.lpot_source.refined_ovsmp_quad_density_discr.ambient_dim + return self.lpot_source.quad_stage2_density_discr.ambient_dim @memoize_method def traversal(self): @@ -390,7 +390,7 @@ class _FMMGeometryData(object): target_info = self.target_info() with cl.CommandQueue(self.cl_context) as queue: - nsources = lpot_src.refined_ovsmp_quad_density_discr.nnodes + nsources = lpot_src.quad_stage2_density_discr.nnodes nparticles = nsources + target_info.ntargets refine_weights = cl.array.zeros(queue, nparticles, dtype=np.int32) @@ -400,7 +400,7 @@ class _FMMGeometryData(object): MAX_LEAF_REFINE_WEIGHT = 32 # noqa tree, _ = code_getter.build_tree(queue, - particles=lpot_src.refined_ovsmp_quad_density_discr.nodes(), + particles=lpot_src.quad_stage2_density_discr.nodes(), targets=target_info.targets, max_leaf_refine_weight=MAX_LEAF_REFINE_WEIGHT, refine_weights=refine_weights, diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index fcb172d2..e3a0575b 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -109,7 +109,7 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): discr_nodes = lpot_source.density_discr.nodes().get(queue) fine_discr_nodes = \ - lpot_source.refined_ovsmp_quad_density_discr.nodes().get(queue) + lpot_source.quad_stage2_density_discr.nodes().get(queue) int_centers = get_centers_on_side(lpot_source, -1) int_centers = np.array([axis.get(queue) for axis in int_centers]) ext_centers = get_centers_on_side(lpot_source, +1) @@ -177,7 +177,7 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): for i, panel_1 in enumerate(iter_elements(lpot_source.density_discr)): for panel_2 in iter_elements(lpot_source.density_discr): check_disk_undisturbed_by_sources(panel_1, panel_2) - for panel_2 in iter_elements(lpot_source.refined_ovsmp_quad_density_discr): + for panel_2 in iter_elements(lpot_source.quad_stage2_density_discr): check_sufficient_quadrature_resolution(panel_1, panel_2) if helmholtz_k is not None: check_panel_size_to_helmholtz_k_ratio(panel_1) -- GitLab