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 ablock
as used by the recursive proxybased skeletonization of the direct solver algorithms. Clusters are represented by aTargetAndSourceClusterList
.
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
.
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]#
 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]#

 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] = 'adaptivelevelrestricted', 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 usingstarts
.
 starts#
An
ndarray
of size(nclusters + 1,)
consisting of nondecreasing integers used to index intoindices
. A cluster \(i\) can be retrieved usingindices[starts[i]:starts[i + 1]]
.
 class pytential.linalg.TargetAndSourceClusterList(targets: IndexList, sources: IndexList)[source]#
Convenience class for working with clusters (subsets) of a matrix.
 nclusters#
 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 alsomake_index_cluster_cartesian_product()
).
 pytential.linalg.make_index_list(indices: ndarray, starts: Optional[ndarray] = None) IndexList [source]#
Wrap a
(indices, starts)
tuple into anIndexList
. 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 twodimensional 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
IntG
s 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 nonconstant 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 selfevaluation.
 class pytential.linalg.direct_solver_symbolic.DOFDescriptorReplacer(source, target)[source]#
Mapper that replaces all the
DOFDescriptor
s 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.