Linear Algebra Routines

Hierarchical Direct Solver

Proxy Point Generation

class pytential.linalg.proxy.ProxyGenerator(source, approx_nproxy=None, ratio=None)[source]
ambient_dim
nproxy

Number of proxy points in a single proxy ball.

source

A pytential.qbx.QBXLayerPotentialSource.

ratio

A ratio used to compute the proxy ball radius. The radius is computed in the \(\ell^2\) norm, resulting in a circle or sphere of proxy points. For QBX, we have two radii of interest for a set of points: the radius \(r_{block}\) of the smallest ball containing all the points and the radius \(r_{qbx}\) of the smallest ball containing all the QBX expansion balls in the block. If the ratio \(\theta \in [0, 1]\), then the radius of the proxy ball is

\[r = (1 - \theta) r_{block} + \theta r_{qbx}.\]

If the ratio \(\theta > 1\), the the radius is simply

\[r = \theta r_{qbx}.\]
ref_points

Reference points on a unit ball. Can be used to construct the points of a proxy ball \(i\) by translating them to center[i] and scaling by radii[i], as obtained by __call__().

__call__(queue, indices, **kwargs)[source]

Generate proxy points for each given range of source points in the discretization in source.

Parameters:
  • queue – a pyopencl.CommandQueue.
  • indices – a sumpy.tools.BlockIndexRanges.
Returns:

a tuple of (proxies, pxyranges, pxycenters, pxyranges), where each element is a pyopencl.array.Array. The sizes of the arrays are as follows: pxycenters is of size (2, nranges), pxyradii is of size (nranges,), pxyranges is of size (nranges + 1,) and proxies is of size (2, nranges * nproxy). The proxy points in a range \(i\) can be obtained by a slice proxies[pxyranges[i]:pxyranges[i + 1]] and are all at a distance pxyradii[i] from the range center pxycenters[i].

pytential.linalg.proxy.partition_by_nodes(discr, use_tree=True, max_nodes_in_box=None)[source]

Generate equally sized ranges of nodes. The partition is created at the lowest level of granularity, i.e. nodes. This results in balanced ranges of points, but will split elements across different ranges.

Parameters:
  • discr – a meshmode.discretization.Discretization.
  • use_tree – if True, node partitions are generated using a boxtree.TreeBuilder, which leads to geometrically close points to belong to the same partition. If False, a simple linear partition is constructed.
  • max_nodes_in_box – passed to boxtree.TreeBuilder.
Returns:

a sumpy.tools.BlockIndexRanges.

pytential.linalg.proxy.partition_by_elements(discr, use_tree=True, max_elements_in_box=None)[source]

Generate equally sized ranges of points. The partition is created at the element level, so that all the nodes belonging to an element belong to the same range. This can result in slightly larger differences in size between the ranges, but can be very useful when the individual partitions need to be resampled, integrated, etc.

Parameters:
  • discr – a meshmode.discretization.Discretization.
  • use_tree – if True, node partitions are generated using a boxtree.TreeBuilder, which leads to geometrically close points to belong to the same partition. If False, a simple linear partition is constructed.
  • max_elements_in_box – passed to boxtree.TreeBuilder.
Returns:

a sumpy.tools.BlockIndexRanges.

pytential.linalg.proxy.partition_from_coarse(resampler, from_indices)[source]

Generate a partition of nodes from an existing partition on a coarser discretization. The new partition is generated based on element refinement relationships in resampler, so the existing partition needs to be created using partition_by_elements(), since we assume that each range contains all the nodes in an element.

The new partition will have the same number of ranges as the old partition. The nodes inside each range in the new partition are all the nodes in resampler.to_discr that were refined from elements in the same range from resampler.from_discr.

Parameters:
  • resampler – a meshmode.discretization.connection.DirectDiscretizationConnection.
  • from_indices – a sumpy.tools.BlockIndexRanges.
Returns:

a sumpy.tools.BlockIndexRanges.

pytential.linalg.proxy.gather_block_neighbor_points(discr, indices, pxycenters, pxyradii, max_nodes_in_box=None)[source]

Generate a set of neighboring points for each range of points in discr. Neighboring points of a range \(i\) are defined as all the points inside the proxy ball \(i\) that do not also belong to the range itself.

Parameters:
  • discr – a meshmode.discretization.Discretization.
  • indices – a sumpy.tools.BlockIndexRanges.
  • pxycenters – an array containing the center of each proxy ball.
  • pxyradii – an array containing the radius of each proxy ball.
Returns:

a sumpy.tools.BlockIndexRanges.

pytential.linalg.proxy.gather_block_interaction_points(source, indices, ratio=None, approx_nproxy=None, max_nodes_in_box=None)[source]

Generate sets of interaction points for each given range of indices in the source discretization. For each input range of indices, the corresponding output range of points is consists of:

  • a set of proxy points (or balls) around the range, which model farfield interactions. These are constructed using ProxyGenerator.
  • a set of neighboring points that are inside the proxy balls, but do not belong to the given range, which model nearby interactions. These are constructed with gather_block_neighbor_points().
Parameters:
Returns:

a tuple (nodes, ranges), where each value is a pyopencl.array.Array. For a range \(i\), we can get the slice using nodes[ranges[i]:ranges[i + 1]].