PoissonSolvers Library Functions

FFTWBasis(domain, ngrid)

A uniform periodic grid of ngrid cells on domain.

xgrid holds ngrid + 1 points, both endpoints included; the last coincides with the first under periodicity and is not a degree of freedom, so length(b) == ngrid.

FFTWDerivative

A derivative of an FFTWSolution, as returned by derivative. Callable, so that derivative(s).(v) broadcasts the way s.(v) does, matching the spline side.

FFTWSolution(basis, coefficients)

The solution on an FFTWBasis or a FiniteDifferenceBasis, evaluated at the nearest grid point.

FiniteDifferenceBasis(domain, ngrid; order = 2)

A uniform periodic grid of ngrid cells on domain, on which $-\Delta$ is the central difference stencil of the given order, 2 or 4.

The grid is the one FFTWBasis lays out: xgrid holds ngrid + 1 points, and the last is not a degree of freedom. This basis selects PoissonSolverMatrixFree.

PoissonSolver(basis)

The solver this basis selects for $-\Delta \phi = \rho$.

PoissonSolverFFT(basis)

A spectral solver for $-\phi'' = \rho$ on a periodic uniform grid.

The transforms and the inverse Laplacian symbol are built once, so solve! allocates nothing when it is given a vector of grid values. Given a function, it allocates the vector of samples it takes first.

The transform scratch is a field of the solver, which is what makes that possible. A solver is therefore not reentrant: two tasks must not call solve! on one solver, even with distinct result vectors. Give each task its own solver.

PoissonSolverMatrixFree(basis::FiniteDifferenceBasis; reltol, maxiter)

A matrix-free solver for $-\phi'' = \rho$ on a periodic uniform grid.

The periodic Laplacian annihilates the constants, so $-\Delta \phi = \rho$ is not well posed. With $R$ the projection onto the constants, $(-\Delta + R) \phi = (1 - R) \rho$ is: the operator is symmetric positive definite, and its solution is the mean-free one, which is the solution PoissonSolverFFT returns. solve! finds it by conjugate gradients, applying the central difference stencil and the projection and building no matrix.

The iteration stops once the residual norm falls to reltol times that of $(1 - R) \rho$, and throws an ErrorException if maxiter iterations do not get it there. reltol bounds the residual, not the error: the error may exceed it by up to the condition number of the operator, which grows as ngrid^2. In Float32 on a fine grid, the solution is therefore accurate to far fewer digits than eps(Float32).

The iteration vectors are fields of the solver, so solve! allocates nothing when it is given a vector of grid values. Given a function, it allocates the vector of samples it takes first. As for PoissonSolverFFT, a solver is therefore not reentrant: give each task its own.

PoissonSolverSpline(basis)

A B-spline Galerkin solver for $-\phi'' = \rho$ on basis.

The stiffness matrix is factorised once, through the representation SimpleSplines chooses for the basis: an FFT for a periodic uniform basis, a banded Cholesky for a Dirichlet one. Both make solve! allocation-free when it is given a coefficient vector. Given a function, it allocates the load vector it samples first — see loadvector.

A periodic stiffness matrix is singular, since $-\phi'' = \rho$ determines the solution only up to a constant. The factorisation deflates that constant, so solve! returns the mean-free solution, and what the operator holds is the assembly itself, which is $O(N)$ sparse.

A periodic factorisation holds transform scratch of its own, which is what makes that possible, so a periodic solver is not reentrant: two tasks must not call solve! on one of them, even with distinct result vectors. The Dirichlet factorisation carries no scratch and does not share this. Give each task its own solver rather than depending on which one it has.

The basis must be periodic or Dirichlet-recombined. Any other basis represents the constants, which leaves the stiffness matrix singular; the constructor rejects it.

Potential(basis)
Potential(basis, rhs)

The solution of $-\Delta \phi = \rho$ on basis, kept together with the solver that produced it so that it can be re-solved in place as rhs changes.

Callable: p(x) evaluates $\phi$, and p(x, d) its d-th derivative. basis(p) and coefficients(p) reach the basis and the coefficient vector; both live inside the solution rather than being stored a second time.

update!(p, rhs) overwrites the coefficients in place. That is visible through p itself because the solution shares its coefficient array rather than copying it — the contract every backend here keeps, and what makes a re-solve from a vector free of allocation.

DirichletBasisSpline(domain, order, ncells)

A B-spline basis of the given order on ncells uniform cells of domain, recombined so that every basis function vanishes at both ends. It spans the solutions of $-\phi'' = \rho$ with homogeneous Dirichlet boundaries, and has two degrees of freedom fewer than the clamped basis.

order is $k = p + 1$, as for PeriodicBasisSpline.

PeriodicBasisSpline(domain, order, ncells)

A periodic B-spline basis of the given order on ncells uniform cells of domain.

order is $k = p + 1$, so order = 2 is the piecewise-linear basis. Note that SimpleSplines constructors take the degree $p$ instead.

PoissonSolution(basis, coefficients)

A callable solution built from basis and coefficients, so that sol(x) evaluates it.

The result shares coefficients rather than copying it, which is what lets update! re-solve in place.

loadvector(p::PoissonSolverSpline, f)

The Galerkin load vector $\int f \phi_i \, dx$ of the function f.

This is not SimpleSplines.l2_projection, which solves with the mass matrix as well. A Poisson right-hand side is the load vector itself.

It returns a fresh vector, so solve! on a function allocates one per call. Build the load vector once and solve with that where the allocation matters.

ndofs(basis)

The number of degrees of freedom of basis.

The backends spell this differently — nbasis for a spline basis, length for a grid — and Potential should not have to know which one it holds.

rhs(p::Potential)

The stored right-hand side buffer of p, which update! solves with by default.

It is a scratch array for the caller to fill, and starts out zero: the right-hand side given to the constructor is solved with, not stored here. So update!(p) on a freshly built p solves for a zero source.

solve!(coefficients, solver, rhs)
solve(solver, rhs)

Solve $-\Delta \phi = \rho$ for the coefficients of $\phi$.

rhs is either the discrete right-hand side or a function, which each backend reduces to one.

solve!(coefficients, solver, rhs)
solve(solver, rhs)

Solve $-\Delta \phi = \rho$ for the coefficients of $\phi$.

rhs is either the discrete right-hand side or a function, which each backend reduces to one.

update!(p::Potential, rhs = rhs(p))

Re-solve p for the right-hand side rhs, in place, and return p.

rhs is either a vector or a function, as for solve!. The new coefficients are written into the array the solution already holds, so the result is visible through p without anything being rebuilt.

With a vector — the stored rhs buffer in particular — nothing is allocated. A function is sampled onto a fresh vector first, which allocates one per call.