PCFFT API

Specifying kernels and points

This package expects points to specified using a struct with field r containing the coordinates of the points. If there are \(N\) points in \(\mathbb{R}^d\) for d=2,3, the r field should be an array with shape [d, N]. It should also contain any other point information, such as a normal vector n with shape [d, N], or a curvature kappa with shape [1,N]. For example, in 3 dimensions the source and target points might be generated by

N = 1000; % number of sources
M = 3000; % number of targets
src_info = struct();
targ_info = struct();
src_info.r = rand(3, N);
src_info.kappa = rand(1, N);
targ_info.r = rand(3, M);
targ_info.n = rand(3, M);

Now, suppose we have a kernel function

\[k : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^{o_1 \times o_2}\]

The PCFFT package requires the kernel to be specified as a function handle or as an object with an eval method. We assume a few things about the kernel function (or its eval method):

The first argument is the source point structure, and the second argument is the target point structure. The kernel function should not require other arguments.

If there are \(N\) source points and \(M\) target points, the kernel function should return an array of shape [opdim(1)*M, opdim(2)*N] containing the kernel evaluations between each target and source point. If the kernel is vector-valued, then the result should be interleaved, i.e. vals(1:opdim(1),1:opdim(2))`` should be the interaction between the first source and the first target.

Here is an example of a kernel function which matches these requirements:

function evals = kern(src_info, targ_info)
    % evaluate the gradient of the electrostatic kernel between N source points src_info.r and
    % M target points targ_info.r, weighted by src_info.kappa
    %
    % Output shape is (3M, N)

    % Shape (M, N)
    rx = src_info.r(1, :) - targ_info.r(1, :).';
    ry = src_info.r(2, :) - targ_info.r(2, :).';
    rz = src_info.r(3, :) - targ_info.r(3, :).';

    dist = sqrt(rx.^2 + ry.^2 + rz.^2);
    gradx = rx ./ dist.^2; gradx = reshape(gradx,1,[]);
    grady = ry ./ dist.^2; grady = reshape(grady,1,[]);
    gradz = rz ./ dist.^2; gradz = reshape(gradz,1,[]);

    % Shape (3M, N)
    evals = reshape([gradx;grady;gradz], 3*size(targ_info.r(:,:),2), size(src_info.r(:,:),2));
    evals = evals .* src_info.kappa(:,:) / (4*pi);
 end

In the following documentation, we will use the types point_info and kernel to refer to structures containing point information and kernel function handles, respectively.

Function API

get_grid(kernel, src_info, targ_info, tol, n_nbr, opts)

Compute the equispaced spreading grid used for the precorrected FFT.

Parameters:

kernel (kernel) – Free-space kernel.
src_info (point_info) – struct giving information about the sources.
targ_info (point_info) – struct describing the target points.
tol (float) – float specifying absolute error tolerance. Error is evaluated at a surface 1.1 * radius of the innermost proxy surface.
n_nbr (int, optional) – int specifying the average number of interactions that must be done directly. Defaults to 1000.
opts (struct, optional) –
options to manipulate the choice of proxy points. Available options:
- opts.multi_shells
  Whether to default to shell-based proxies. Defaults to false. Accelerates the precomputation for kernels that are known to not satisfy Green’s identity
- opts.proxy_der
  Number of radial derivatives to use in the proxy Can be a number between 0 and 2. Defaults to 0. This option can avoid invoking shells for kernels that are derived from high order PDEs. (see wrap_kern_der)
- opts.halfside
  Manually set box size (and implicitly n_nbr). Only recommended for expert users. (See pcff_fmm3dbie_demo.m) Useful for plotting BIE solutions without overly small boxes.

Returns:

grid_info (GridInfo) – GridInfo object specifying the regular grid used for spreading.
proxy_info (ProxyInfo) – ProxyInfo object specifying the proxy surface(s) used in the proxy point method.

get_spread(kern_0, kern_der, src_info, grid_info, proxy_info, der_fields)

This routine returns the matrix that maps charge strengths at srcinfo.r to charge strengths on the equispaced grid.

Parameters:

kern_0 (kernel) – The free-space kernel, which must be scalar-valued
kern_der (kernel) – Some derivative of the free-space kernel. This can be left empty by passing in an empty array, in which case the free-space kernel will be used. Each pairwise interaction must of shape [1, opdim]
src_info (point_info) – Specifies the source points.
grid_info (GridInfo) – object describing the regular grid
proxy_info (ProxyInfo) – object describing the proxy points
der_fields (cell array, optional) – Contains field names that must be attached to the source point in kern_der, e.g. {‘r’, ‘n’, ‘kappa’}. If omitted or left empty this argument defaults to {‘r’}.

Returns:

A_spread (sparse matrix [nreg, opdim*nsrc]) – Maps source strengths to equivalent strengths on the regular grid.
sort_info (SortInfo) – Object describing the sorting of source points into bins

get_addsub(kern_0, kern_st, grid_info, proxy_info, sort_info_s, sort_info_t, A_spread_s, A_spread_t)

Compute the correction for near-field interactions.

Parameters:

kern_0 (kernel) – Free-space kernel, which must be scalar-valued
kern_st (kernel) – Direct interaction kernel, which must be a combination of derivatives of the free-space kernel. If left empty the free-space kernel will be used. Each pairwise interaction must of shape [opdim(1), opdim(2)]
grid_info (GridInfo) – GridInfo object describing the regular grid.
proxy_info (ProxyInfo) – ProxyInfo object describing the proxy points.
sort_info_s (SortInfo) – Specifies how source points are sorted into bins.
sort_info_t (SortInfo) – Specifies how target points are sorted into bins.
A_spread_s (sparse matrix [nreg, opdim(2)*nsrc]) – Maps source strengths to equivalent strengths on the regular grid.
A_spread_t (sparse matrix [nreg, opdim(1)*ntarg]) – Maps target strengths to equivalent strengths on the regular grid. Its adjoint is used to map regular grid strengths to equivalent strengths on the target points.

Returns:

A_addsub (sparse matrix [n_targ, n_src]) – A sparse matrix which corrects for the incorrect near-field interactions computed using the spreading matrices.

get_kernhat(kern_0, grid_info)

Evaluate Fourier transform of free-space kernel on the grid in grid_info. Ignores the self-interaction.

Parameters:

kern_0 (kernel) – Free-space kernel.
grid_info (GridInfo) – object describing the regular grid.

Returns:

kern_hat (matrix [2*ngrid(dim), …, 2*ngrid(1)]) – Fourier transform of the kernel evaluated on the grid.

pcfft_apply(mu, A_spread_s, A_spread_t, A_addsub, kern_0hat)

Compute N-body sum using a precorrected FFT.

Parameters:

mu (matrix [opdim(2)*nsrc, 1]) – Source strengths.
A_spread_s (sparse matrix) – source spreading matrix (see get_spread).
A_spread_t (sparse matrix) – target spreading matrix.
A_addsub (sparse matrix) – near field corrections matrix (see get_addsub).
kern_0hat (matrix) – Fourier transform of background kernel (see get_kernhat).

Returns:

u (matrix) – potential.

Built-in classes

class GridInfo

Describes the regular grid used for spreading; is created by get_grid().

Stores information about the regular grid used by the spreading routines. The grid is divided into non-overlapping spreading bins, each of which is spread to a larger spreading box independently.

Variables:

ngrid (array [dim, 1]) – number of regular grid points per dimension.
Lbd (array [dim, 2]) – containing the points (xmin, ymin, zmin) and (xmax, ymax, zmax). Describes the bounding box of the union of source and target points.
dx (float) – grid spacing.
nspread (int) – number of spreading points per dimension required for the requested accuracy.
nbinpts (int) – width in dx’s of the spreading bin.
rpad (int) – Width of the margin between the bounding box specified by Lbd and the regular grid points.
r (array [dim, npts]) – padded regular grid points.
dim (int) – problem dimension.
nbin (array [dim, 1]) – number of spreading bins per dimension.
offset (int) – Number of dx’s used to pad the grid.
rmax (array [dim, 1]) – maximum coordinate value of the padded grid.
rmin (array [dim, 1]) – minimum coordinate value of the padded grid.
n_nbr (int) – average number of near-field neighbours.
zero_bin (array [dim, nbinpts^dim]) – Grid points of a spreading bin centered at the origin.
zero_box (array [dim, nspread^dim]) – Grid points of a spreading box centered at the origin.
center_bin (int) –

Linear index of a bin approximately at the center of the grid.
Computed as floor(nbin(1)/2) * nbin(2) + floor(nbin(2)/2).

nbr_offsetsarray [n_nbr, 1]
Used for indexing neighboring bins. n_nbr is the number of spreading bins which are treated directly.
idx_cutoff (int) – The number of bins in each direction which are treated directly.

class ProxyInfo

Describes the proxy surface for spreading weights.

Stores metadata describing the proxy surface used by the spreading routines to compute spreading weights. The proxy surface can consist of one or more concentric shells of points.

Variables:

dim (int) – problem dimension (integer).
n_points_total (int) – total number of proxy points = nproxy * nshell.
nproxy (int) – number of proxy points per shell.
nshell (int) – number of concentric shells in the proxy.
halfside (float) – Specifies the halfside of the spreading box.
crad (float) – Parameter for determining the radius of the innermost proxy shell.
tol (float) – relative L_inf error tolerance for generating the proxy surface.
radius (float) – Radius of the innermost proxy shell. Used for determining near-field neighbors. = sqrt(dim) * halfside * crad.
r (array [dim, n_points_total]) – array of proxy point coords for a spreading box centered at the origin.
proxy_der (int) – number of derivatives used in proxy compression. See wrap_kern_der.

class SortInfo

Information about sorting of source points into bins.

The source points are sorted into bins, and various sorted arrays and indices are stored to avoid repeated binsort operations.

Variables:

r_srt (array [dim, n_src]) – Contains the source points sorted by bin index.
binid_srt (array [n_src, 1]) – Contains the bin indexes of the sorted source points.
ptid_srt (array [n_src, 1]) – Contains the indices of the sorting of the source points. So r_srt = src_info.r(:, ptid_srt).
id_start (array [N_bins + 1, 1]) – Contains the indices of the start of each bin in r_srt.
data_srt (struct) – struct of data associated with each point in r_srt. Each field must be of shape [l, n_src] for some l.