Linear Algebra Routines#

In the linear algebra parts of pytential, the following naming scheme is used:

block refers to a piece of a vector operator, e.g. the \(S_{xx}\) component of the Stokeslet.
cluster refers to a piece of a block as used by the recursive proxy-based skeletonization of the direct solver algorithms. Clusters are represented by a TargetAndSourceClusterList.

GMRES#

pytential.linalg.gmres.gmres(op: Callable[[ArrayOrContainerT], ArrayOrContainerT], rhs: ArrayOrContainerT, restart: Optional[int] = None, tol: Optional[float] = None, x0: Optional[ArrayOrContainerT] = None, inner_product: Optional[Callable[[ArrayOrContainerT, ArrayOrContainerT], float]] = None, maxiter: Optional[int] = None, hard_failure: Optional[bool] = None, no_progress_factor: Optional[float] = None, stall_iterations: Optional[int] = None, callback: Optional[Callable[[ArrayOrContainerT], None]] = None, progress: bool = False, require_monotonicity: bool = True) → GMRESResult[source]#

Solve a linear system \(Ax = b\) using GMRES with restarts.

Parameters:

op – a callable to evaluate \(A(x)\).
rhs – the right hand side \(b\).
restart – the maximum number of iteration after which GMRES algorithm needs to be restarted
tol – the required decrease in residual norm (relative to the rhs).
x0 – an initial guess for the iteration (a zero array is used by default).
inner_product – a callable with an interface compatible with numpy.vdot() that returns a host scalar.
maxiter – the maximum number of iterations permitted.
hard_failure – if True, raise GMRESError in case of failure.
stall_iterations – number of iterations with residual decrease below no_progress_factor indicates stall. Set to 0 to disable stall detection.

Returns:

a GMRESResult.

class pytential.linalg.gmres.GMRESResult(solution: ArrayContainer, residual_norms: Sequence[float], iteration_count: int, success: bool, state: str)[source]#

solution#

residual_norms#

iteration_count#

success#: A bool indicating whether the iteration succeeded.

state#: A description of the outcome.

exception pytential.linalg.gmres.GMRESError[source]#

class pytential.linalg.gmres.ResidualPrinter(inner_product=<function structured_vdot>)[source]#

Hierarchical Direct Solver#

Warning

All the classes and routines in this module are experimental and the API can change at any point.

Proxy Point Generation#

class pytential.linalg.ProxyPointSource(lpot_source: QBXLayerPotentialSource, proxies: ndarray)[source]#

get_expansion_for_qbx_direct_eval(base_kernel, target_kernels)[source]#: Wrapper around pytential.qbx.QBXLayerPotentialSource.get_expansion_for_qbx_direct_eval to allow this class to be used by the matrix builders.

class pytential.linalg.ProxyPointTarget(lpot_source: QBXLayerPotentialSource, proxies: ndarray)[source]#

class pytential.linalg.ProxyClusterGeometryData(places: GeometryCollection, dofdesc: DOFDescriptor, srcindex: IndexList, pxyindex: IndexList, points: ndarray, centers: ndarray, radii: ndarray, _cluster_radii: Optional[ndarray] = None)[source]#

srcindex#: A IndexList describing which cluster of points each proxy ball was created from.

pxyindex#: A IndexList describing which proxies belong to which cluster.

points#: A concatenated array of all the proxy points. Can be sliced into using pxyindex (shape (dim, nproxies)).

centers#: An array of all the proxy ball centers (shape (dim, nclusters)).

radii#: An array of all the proxy ball radii (shape (nclusters,)).

nclusters#

__init__(places: GeometryCollection, dofdesc: DOFDescriptor, srcindex: IndexList, pxyindex: IndexList, points: ndarray, centers: ndarray, radii: ndarray, _cluster_radii: Optional[ndarray] = None) → None#

to_numpy(actx: PyOpenCLArrayContext) → ProxyClusterGeometryData[source]#

as_sources() → ProxyPointSource[source]#

as_targets() → ProxyPointTarget[source]#

class pytential.linalg.ProxyGeneratorBase(places: GeometryCollection, approx_nproxy: Optional[int] = None, radius_factor: Optional[float] = None, norm_type: str = 'linf', _generate_ref_proxies: Optional[Callable[[int], ndarray]] = None)[source]#

nproxy#: Number of proxy points in a single proxy ball.

__init__(places: GeometryCollection, approx_nproxy: Optional[int] = None, radius_factor: Optional[float] = None, norm_type: str = 'linf', _generate_ref_proxies: Optional[Callable[[int], ndarray]] = None) → None[source]#

Parameters:

approx_nproxy – desired number of proxy points. In higher dimensions, it is not always possible to construct a proxy ball with exactly this number of proxy points. The exact number of proxy points can be retrieved with nproxy.
radius_factor – Factor multiplying the cluster radius (i.e radius of the bounding box) to get the proxy ball radius.
norm_type – type of the norm used to compute the centers of each cluster. Supported values are "linf" and "l2".

__call__(actx: PyOpenCLArrayContext, source_dd: Optional[pytential.symbolic.dof_desc.DOFDescriptorLike], dof_index: IndexList, **kwargs: Any) → ProxyClusterGeometryData[source]#

Generate proxy points for each cluster in dof_index_set with nodes in the discretization source_dd.

Parameters:: source_dd – a DOFDescriptor for the discretization on which the proxy points are to be generated.

class pytential.linalg.ProxyGenerator(places: GeometryCollection, approx_nproxy: Optional[int] = None, radius_factor: Optional[float] = None, norm_type: str = 'linf', _generate_ref_proxies: Optional[Callable[[int], ndarray]] = None)[source]#

A proxy point generator that only considers the points in the current cluster when determining the radius of the proxy ball.

Inherits from ProxyGeneratorBase.

class pytential.linalg.QBXProxyGenerator(places: GeometryCollection, approx_nproxy: Optional[int] = None, radius_factor: Optional[float] = None, norm_type: str = 'linf', _generate_ref_proxies: Optional[Callable[[int], ndarray]] = None)[source]#

A proxy point generator that also considers the QBX expansion when determining the radius of the proxy ball.

Inherits from ProxyGeneratorBase.

pytential.linalg.partition_by_nodes(actx: PyOpenCLArrayContext, places: GeometryCollection, *, dofdesc: Optional[pytential.symbolic.dof_desc.DOFDescriptorLike] = None, tree_kind: Optional[str] = 'adaptive-level-restricted', max_particles_in_box: Optional[int] = None) → IndexList[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:

dofdesc – a DOFDescriptor for the geometry in places which should be partitioned.
tree_kind – if not None, it is passed to boxtree.TreeBuilder.
max_particles_in_box – value used to control the number of points in each partition (and thus the number of partitions). See the documentation in boxtree.TreeBuilder.

pytential.linalg.gather_cluster_neighbor_points(actx: PyOpenCLArrayContext, pxy: ProxyClusterGeometryData, *, max_particles_in_box: Optional[int] = None) → IndexList[source]#: Generate a set of neighboring points for each cluster of points in pxy. Neighboring points of a cluster \(i\) are defined as all the points inside the proxy ball \(i\) that do not also belong to the cluster itself.

Misc#

class pytential.linalg.IndexList(indices: ndarray, starts: ndarray)[source]#

Convenience class for working with clusters (subsets) of an array.

nclusters#

indices#: An ndarray of not necessarily continuous or increasing integers representing the indices of a global array. The individual cluster slices are delimited using starts.

starts#: An ndarray of size (nclusters + 1,) consisting of nondecreasing integers used to index into indices. A cluster \(i\) can be retrieved using indices[starts[i]:starts[i + 1]].

cluster_size(i: int) → int[source]#

cluster_indices(i: int) → ndarray[source]#

Returns:: a view into the indices array for the range corresponding to cluster i.

cluster_take(x: ndarray, i: int) → ndarray[source]#

Returns:: a subset of x corresponding to the indices in cluster i. The returned array is a copy (not a view) of the elements of x.

class pytential.linalg.TargetAndSourceClusterList(targets: IndexList, sources: IndexList)[source]#

Convenience class for working with clusters (subsets) of a matrix.

nclusters#

targets#: An IndexList encapsulating target cluster indices.

sources#: An IndexList encapsulating source cluster indices.

cluster_shape(i: int, j: int) → Tuple[int, int][source]#

Returns:: the shape of the cluster (i, j), where i indexes into the targets and j into the sources.

cluster_indices(i: int, j: int) → Tuple[ndarray, ndarray][source]#

Returns:: a view into the indices that make up the cluster (i, j).

cluster_take(x: ndarray, i: int, j: int) → ndarray[source]#

Returns:: a subset of the matrix x corresponding to the indices in the cluster (i, j). The returned array is a copy of the elements of x.

flat_cluster_take(x: ndarray, i: int) → ndarray[source]#

Returns:: a subset of an array x corresponding to the indices in the cluster i. Unlike cluster_take(), this method indexes into flattened cluster arrays (see also make_index_cluster_cartesian_product()).

pytential.linalg.make_index_list(indices: ndarray, starts: Optional[ndarray] = None) → IndexList[source]#

Wrap a (indices, starts) tuple into an IndexList.

Parameters:: starts – if None, then indices is expected to be an object array of indices, so that the starts can be reconstructed.

pytential.linalg.make_index_cluster_cartesian_product(actx: PyOpenCLArrayContext, mindex: TargetAndSourceClusterList) → Tuple[Array, Array][source]#

Constructs a cluster by cluster Cartesian product of all the indices in mindex.

The indices in the resulting arrays are laid out in C order. Retrieving two-dimensional data for a cluster \(i\) using the resulting index arrays can be done as follows

offsets = np.cumsum([0] + [
    mindex.targets.cluster_size(i) * mindex.sources.cluster_size(i)
    for i in range(mindex.nclusters)
    ])
istart = offsets[i]
iend = offsets[i + 1]

cluster_1d = x[tgtindices[istart:iend], srcindices[istart:iend]]
cluster_2d = cluster_1d.reshape(*mindex.cluster_shape(i, i))

assert np.allclose(cluster_2d, mindex.cluster_take(x, i, i))

The result is equivalent to cluster_take(), which takes the Cartesian product as well.

Returns:: a tuple containing (tgtindices, srcindices), where the type of the arrays is the base array type of actx.

Internal Functionality#

class pytential.linalg.direct_solver_symbolic.KernelTransformationRemover[source]#

A mapper that removes the transformations from the kernel of all IntGs in the expression.

This includes source and target derivatives and other such transformations. Any unnecessary kernel arguments are also removed from kernel_arguments.

This mapper is meant to be used in the directs solver for proxy interaction, where it is not possible to evaluate source or target directional derivatives.

class pytential.linalg.direct_solver_symbolic.IntGTermCollector[source]#

A mapper that removes all non-IntG terms from the expression and all their non-constant factors.

In particular, an expression of the type

\[\sum_{i = 0}^N f_i(\mathbf{x}, \sigma) + \sum_{i = 0}^M c_i g_i(\mathbf{x}) \mathrm{IntG}_i(\mathbf{x})\]

is reduced to

\[\sum_{i = 0}^M c_i \mathrm{IntG}_i(\mathbf{x}).\]

The intended used of this transformation is to allow the evaluation of the proxy interactions in the direct solver for a given expression meant for self-evaluation.

class pytential.linalg.direct_solver_symbolic.DOFDescriptorReplacer(source, target)[source]#

Mapper that replaces all the DOFDescriptors in the expression with the given ones.

This mapper is meant to allow for evaluation of proxy interactions in the direct solver when the given expression is already partially (or fully) tagged.

__init__(source, target)[source]#

Parameters:

source – a descriptor for all expressions to be evaluated on the source geometry.
target – a descriptor for all expressions to be evaluate on the target geometry.