diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 50fb8e5887266235ac7864fee91986217e06be96..f0a0a6fb2ad1f2a57451088eb7cafcd00d1f8ce9 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,39 +238,70 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): # }}} @property - def base_fine_density_discr(self): + def stage2_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_interp_to_ovsmp_quad_connection(self): + from meshmode.discretization.connection import make_same_mesh_connection + + return make_same_mesh_connection( + self.quad_stage2_density_discr, + self.stage2_density_discr) + + @property + @memoize_method + def quad_stage2_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.stage2_density_discr.mesh, 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: + # 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.stage2_density_discr, + p.area_element(self.ambient_dim, self.dim) + )(queue)) + + qweight = bind(self.quad_stage2_density_discr, p.QWeight())(queue) + + return (area_element.with_queue(queue)*qweight).with_queue(None) + + # }}} + @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.fine_density_discr, self.base_fine_density_discr) + conn = self.refined_interp_to_ovsmp_quad_connection - 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 +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.base_fine_density_discr) + return utils.element_centers_of_mass(self.stage2_density_discr) @memoize_method def _expansion_radii(self, last_dim_length): @@ -364,7 +396,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.stage2_density_discr, last_dim_length) @memoize_method def qbx_fmm_geometry_data(self, target_discrs_and_qbx_sides): @@ -646,7 +678,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): evt, output_for_each_kernel = lpot_applier( queue, target_discr.nodes(), - self.fine_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])) @@ -659,7 +691,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.quad_stage2_density_discr.nodes(), [strengths], **kernel_args) qbx_forced_limit = o.qbx_forced_limit @@ -708,7 +741,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): lpot_applier_on_tgt_subset( queue, targets=target_discr.nodes(), - sources=self.fine_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 96bacdf18deeadcdd0db9e2131883a52ccff69a4..9a2914ae28d5db325232166ffe3a87f7ea8d61aa 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.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.fine_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.fine_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.fine_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 6c47a1047e735bfa27d960129d96930dddde852a..e1b745634104d5670753257ab998728a7eb0fa9e 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_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_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.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 = [] - base_fine_density_discr = lpot_source.density_discr + stage2_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_stage2_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_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( - base_fine_density_discr, + stage2_density_discr, refiner, refine_flags, group_factory, debug) - base_fine_density_discr = conn.to_discr + stage2_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 edf92dfe1af8a1af70042b761999eed0ed8a6eb2..c27b9a870e2d83e1fee45282bb5b7fd4d8a2bf47 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_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_base_fine_discr=use_base_fine_discr) + use_stage2_discr=use_stage2_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_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.base_fine_density_discr`` + ``lpot_source.density_discr`` or + ``lpot_source.stage2_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_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: @@ -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_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_base_fine_discr + if not use_stage2_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_stage2_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_stage2_discr: + density_discr = lpot_source.quad_stage2_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_stage2_discr else None) # Transfer all tree attributes. diff --git a/pytential/source.py b/pytential/source.py index 75d0db70ee6babe762f9589f67fc8da877fd2379..f6cabf2087859686ce1a6975c51c6c26d3e8d267 100644 --- a/pytential/source.py +++ b/pytential/source.py @@ -151,7 +151,8 @@ class LayerPotentialSourceBase(PotentialSource): .. rubric:: Discretizations .. attribute:: density_discr - .. attribute:: fine_density_discr + .. attribute:: stage2_density_discr + .. attribute:: quad_stage2_density_discr .. attribute:: resampler .. method:: with_refinement @@ -260,27 +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: - # fine_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.fine_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 2012b9928c0dad7e10947ba110c5ca65e9041faa..9f2f101ae64c8a654ec28b41b98fc2b92f395d3b 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.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 fcd5298e798c41bc9853c2d64491d0e84fce2c3c..43ad6e59198925e76a80306a0af532fd382d53cc 100644 --- a/pytential/unregularized.py +++ b/pytential/unregularized.py @@ -88,8 +88,28 @@ 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: + # 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.resampler( + queue, + bind( + self.density_discr, + p.area_element(self.ambient_dim, self.dim) + )(queue)) + + qweight = bind(self.quad_stage2_density_discr, p.QWeight())(queue) + + return (area_element.with_queue(queue)*qweight).with_queue(None) + @property - def fine_density_discr(self): + def quad_stage2_density_discr(self): return self.density_discr def resampler(self, queue, f): @@ -160,7 +180,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.quad_stage2_density_discr.nodes(), [strengths], **kernel_args) result.append((o.name, output_for_each_kernel[o.kernel_index])) @@ -342,11 +363,11 @@ class _FMMGeometryData(object): @property def coord_dtype(self): - return self.lpot_source.fine_density_discr.nodes().dtype + return self.lpot_source.quad_stage2_density_discr.nodes().dtype @property def ambient_dim(self): - return self.lpot_source.fine_density_discr.ambient_dim + return self.lpot_source.quad_stage2_density_discr.ambient_dim @memoize_method def traversal(self): @@ -369,7 +390,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.quad_stage2_density_discr.nnodes nparticles = nsources + target_info.ntargets refine_weights = cl.array.zeros(queue, nparticles, dtype=np.int32) @@ -379,7 +400,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.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 12c816f4a15d771a1ff94502e4bb516ba5d2999c..e3a0575b8fea84a468ac580c5fb0b51d38ede377 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.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) @@ -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.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)