Particle FMM in C++

This tutorial walks through examples/src/example1.cpp, which evaluates the gradient of the Laplace potential of N random particle sources at N random targets, then verifies the result against a direct \(O(N^2)\) sum.

Build and run it with:

make examples/bin/example1
mpirun -n 2 ./examples/bin/example1 -N 100000 -m 10 -omp 4

1. Select a kernel

Every kernel is available as a pvfmm::Kernel<Real> descriptor (Kernel functions); ker_dim[0]/ker_dim[1] give the source/target dimensions per point:

const pvfmm::Kernel<double>& kernel_fn = pvfmm::LaplaceKernel<double>::gradient();

2. Set up sources and targets

Coordinates live in \([0,1]^3\) in array-of-structures layout. The example draws uniform random points (point_distrib is an example helper) and random densities; it uses single-layer sources only (the double-layer arrays are created empty but could carry dipole sources — see C++ interface):

vec trg_coord = point_distrib<double>(RandUnif, N, comm);
vec  sl_coord = point_distrib<double>(RandUnif, N, comm);
vec  dl_coord = point_distrib<double>(RandUnif, 0, comm);   // no double-layer sources

vec sl_den(n_sl * kernel_fn.ker_dim[0]);                    // 1 value/pt for Laplace
vec dl_den(n_dl * (kernel_fn.ker_dim[0] + PVFMM_COORD_DIM)); // density + normal
for (auto& a : sl_den) a = drand48();

3. Build the tree and load the operators

PtFMM_CreateTree partitions the points across MPI ranks and builds the distributed octree (max_pts bounds the number of source points per leaf). PtFMM::Initialize loads/precomputes the translation operators for the kernel at the requested multipole order, and SetupFMM binds them to the tree:

size_t max_pts = 600;
auto* tree = PtFMM_CreateTree(sl_coord, sl_den, dl_coord, dl_den,
                              trg_coord, comm, max_pts, pvfmm::FreeSpace);

pvfmm::PtFMM<double> matrices;
matrices.Initialize(mult_order, comm, &kernel_fn);

tree->SetupFMM(&matrices);

4. Evaluate

vec trg_value;
PtFMM_Evaluate(tree, trg_value, n_trg);

trg_value receives n_trg * ker_dim[1] values (3 per target for the Laplace gradient), partitioned across ranks according to n_trg.

5. Re-evaluate with new densities

When only the densities change (same particle positions), skip the tree and setup work — clear the FMM data and pass the new densities directly:

tree->ClearFMMData();
// ... modify sl_den / dl_den ...
PtFMM_Evaluate(tree, trg_value, n_trg, &sl_den, &dl_den);

6. Verify and clean up

The example samples a subset of targets, computes the direct sum with the same kernel (kernel_fn.ker_poten) gathered over all ranks, and prints the maximum absolute/relative error. The error converges exponentially in the multipole order — measured for this example: 1e-3 (m=4), 1e-5 (m=6), 3e-7 (m=8), 6e-9 (m=10), 1e-10 (m=12) — about two digits per increment of 2 in m. Finally:

delete tree;

The same workflow from the other language interfaces: Using the C and Fortran interfaces, Using the Python and Julia bindings.