Galerkin projection of delta functions

[michakraus]

We are using BSplineKit in particle codes and repeatedly require to project particle distributions onto a spline basis. Effectively that amounts to a L2 projection of sums of delta functions, i.e., if a spline is given by

$$ \phi (x) = \sum \limits_{i=1}^{n} \varphi_i \psi_i (x) $$

we have the following L2 projection

$$ \int \psi{i} (x) \phi (x) dx = \int \psi{i} (x) \sum \limits{a=1}^{N} \delta ( x - x{a} ) dx \quad \forall i $$

and thus the degrees of freedom are given by

$$ \varphii = M{ij}^{-1} \sum \limits_{a=1}^{N} \int \psij (x{a}) \quad \forall i $$

What would be a good and efficient way to implement this in BSplineKit?

Currently we evaluate the right-hand side using evaluate_all, but due to the required index arithmetics when distributing the contribution of one particle to the degrees of freedom of the various basis functions, in whose support the particle resides, this is not particularly nice, especially with periodic splines.

Is there a better way to implement this? We'd be happy to contribute and add this functionality to BSplineKit if this makes sense, but we may need a little guidance.

[jipolanco]

Hi, I'm guessing your last equation should be something like:

$$ M_{ij} \varphij = \sum{\alpha = 1}^N \psii(x\alpha) $$

In that case yes, evaluate_all would be the right way to evaluate the right-hand side. I think it would be interesting to have a convenient way of computing this projection and distributing the contribution of each $x_\alpha$.

I will try to come up with an example of how to do this right now for periodic bases, which should help define a high-level interface for this.

[michakraus]

Yes, of course. I forgot the inverse on the mass matrix. Just added that.

An example would be great. We should be able to work our way from there.

[jipolanco]

This seems to work:

using BSplineKit

N = 100
L = 1
breaks = range(0, 1; length = N + 1)[1:N]
B = PeriodicBSplineBasis(BSplineOrder(4), breaks, L)

# Particle locations
xs = rand(5)

# Compute right-hand side
fs = Splines.PeriodicVector(zeros(length(B)))
for x ∈ xs
    ilast, bs = B(x)  # same as `evaluate_all`
    # Iterate over evaluated basis functions.
    # The indices of the evaluated basis functions are ilast:-1:(ilast - k + 1),
    # where k is the spline order.
    for (δi, bi) ∈ pairs(bs)
        i = ilast + 1 - δi
        fs[i] += bi
    end
end

M = galerkin_matrix(B)
cs = M \ parent(fs)  # get coefficients
S = Spline(B, cs)    # construct spline from coefficients

## Plots

using CairoMakie

fig = Figure()
ax = Axis(fig[1, 1]; xgridvisible = false, ygridvisible = false)
lines!(ax, 0..1, x -> S(x))
vlines!(ax, xs; color = :grey)
save("test.svg", fig)

Note that in the specific case of periodic splines, it's more convenient to work with the internal PeriodicVector type, which just wraps a regular vector allowing for periodic (cyclic) indexing.

[jipolanco]

Just for fun, note that you can also compute an approximation of the RHS using integration by parts and the galerkin_projection function:

fs_alt = galerkin_projection(B, Derivative(1)) do x
    count(y -> x ≥ y, xs)  # integral of the sum of deltas
end
cs_alt = -1 .* (M \ fs_alt)
S_alt = Spline(B, cs_alt)

This is just an approximation because galerkin_projection uses Gauss-Legendre quadratures, which don't capture the discontinuity in the integral of each delta. But it works quite well:

fig = Figure()
ax = Axis(fig[1, 1]; xgridvisible = false, ygridvisible = false)
lines!(ax, 0..1, x -> S(x))
lines!(ax, 0..1, x -> S_alt(x); linestyle = :dash)
vlines!(ax, xs; color = :grey)
save("test.svg", fig)

[michakraus]

Many thanks for the detailed example, @jipolanco ! This was very helpful.

We are integrating this now in our packages PoissonSolvers.jl and ParticleMethods.jl.