Gravitational Example

FastMultipole is designed to incorporate easily into your existing Julia code with minimal effort. In this section and the Vortex Filament Example, we demonstrate how this can be done.

First, we'll review the interface functions used by the gravitational point mass model used in Quick Start. This code can also be found under FastMultipole/test/gravitational.jl.

To better understand how the FastMultipole interface functions, let's take a look at the data structures we'll use to define our point masses:

using FastMultipole
using FastMultipole.StaticArrays

const i_POTENTIAL = 1:1   # index of the gravitational potential
const i_GRADIENT = 5:7    # index of the gravitational acceleration
const i_HESSIAN = 8:16    # index of the hessian matrix

# a single point mass
struct Body{TF}
    position::SVector{3,TF}
    radius::TF
    strength::TF
end

# container for a system of `Body`'s
struct Gravitational{TF}
    bodies::Vector{Body{TF}}
    potential::Matrix{TF}
end

# constructor
function Gravitational(bodies::Matrix)
    nbodies = size(bodies)[2]
    bodies2 = [Body(SVector{3}(bodies[1:3,i]),bodies[4,i],bodies[5,i]) for i in 1:nbodies]
    potential = zeros(eltype(bodies), 16, nbodies)
    return Gravitational(bodies2,potential)
end

Main.Gravitational

In short, Body objects represent individual point masses, and the Gravitational object contains a group of them, along with a matrix of the potential, velocity, and velocity gradient at each body. The observant reader will notice that the strength field of Body represents its mass.

Overloading body_to_multipole!

The function body_to_multipole! is used to generate multipole expansions for our particular system. This is done be overloading FastMultipole.body_to_multipole! for our data structure. Convenience functions exist within FastMultipole to simplify this complicated function into a one-liner:

FastMultipole.body_to_multipole!(system::Gravitational, args...) =
    FastMultipole.body_to_multipole!(Point{Source}, system, args...)

The ::Gravitational type annotation is the key that allows FastMultipole to dispatch on ::Gravitational systems. Point{Source} is used to indicate that our bodies are points and induce a source (i.e. $1/r$) potential. Other convenience functions exist in FastMultipole for any combination of Point, Filament, or Panel geometries using Source, Dipole, or Vortex kernels, and are specified when overloading body_to_multipole! as <geometry>{<kernel>}. For example, if my model used constant vortex sheet panels, I would replace Point{Source} in the example above to Panel{Vortex}.

We use the fast recursive method of generating exact coefficients for panels as developed by [1]. We have derived our own formulae for vortex filaments and panels, which are in the process of publication.

FastMultipole also requires the user to indicate if a source system induces a vector potential. This is done by overloading the has_vector_potential function as follows:

function FastMultipole.has_vector_potential(system::Gravitational)
    return false
end

This will allow FastMultipole to avoid unnecessary calculations if the vector potential is not needed. Since Gravitational systems do not induce a vector potential, we return false.

Buffer Interface Functions

FastMultipole allocates a buffer matrix in which to sort and store key system information. The buffer interface is created by overloading several functions for the user-defined system type. The first of these functions is FastMultipole.source_system_to_buffer!, which populates a matrix with the important information about each body, column-wise. Overloading for ::Gravitational systems looks like this:

function FastMultipole.source_system_to_buffer!(buffer, i_buffer, system::Gravitational, i_body)
    buffer[1:3, i_buffer] .= system.bodies[i_body].position
    buffer[4, i_buffer] = system.bodies[i_body].radius
    buffer[5, i_buffer] = system.bodies[i_body].strength
end

Here, the i_bodyth body position, radius, and strength are copied to the correct positions in the i_bufferth column of the buffer matrix.

It is worth noting that there are ways of defining these functions that would harm performance, e.g. by allocating an array each time the velocity is requested. If you notice FastMultipole's performance is poor, this and the other interface functions are good places to check.

Because every system can have a unique buffer structure, we need to overload a few additional functions to tell FastMultipole how to allocate and access it. The first function, data_per_body, returns the number of rows in the buffer matrix that are used to define a single body. The second function, get_position, returns the position of the ith body in the system. The third function, strength_dims, returns the number of rows in the buffer matrix that are used to define the strength of a single body. Finally, we overload get_n_bodies to return the number of bodies in our system. For our Gravitational system, we define:

Base.eltype(::Gravitational{TF}) where TF = TF

function FastMultipole.data_per_body(system::Gravitational)
    return 5
end

function FastMultipole.get_position(system::Gravitational, i)
    return system.bodies[i].position
end

function FastMultipole.strength_dims(system::Gravitational)
    return 1
end

FastMultipole.get_n_bodies(system::Gravitational) = length(system.bodies)

These functions provide enough information for FastMultipole to allocate the buffer matrix and access the data it needs.

Finally, we need to tell FastMultipole how to transfer the results of the FMM call back to the user-defined system. This is done by overloading the FastMultipole.buffer_to_target_system! function. Overloading for ::Gravitational systems looks like this:

function FastMultipole.buffer_to_target_system!(target_system::Gravitational, i_target, ::FastMultipole.DerivativesSwitch{PS,VS,GS}, target_buffer, i_buffer) where {PS,VS,GS}
    # get values
    TF = eltype(target_buffer)
    scalar_potential = PS ? FastMultipole.get_scalar_potential(target_buffer, i_buffer) : zero(TF)
    velocity = VS ? FastMultipole.get_gradient(target_buffer, i_buffer) : zero(SVector{3,TF})
    hessian = GS ? FastMultipole.get_hessian(target_buffer, i_buffer) : zero(SMatrix{3,3,TF,9})

    # update system
    target_system.potential[i_POTENTIAL[1], i_target] = scalar_potential
    target_system.potential[i_GRADIENT, i_target] .= velocity
    for (jj,j) in enumerate(i_HESSIAN)
        target_system.potential[j, i_target] = hessian[jj]
    end
end

We note that the convenience functions get_gradient and get_hessian return ::SVector{3} and ::SMatrix{3,3}, respectively, to reduce allocations.

Overloading direct!

The last required interface function is FastMultipole.direct!, which evaluates the potential at a target system using a source system without multipole acceleration. This is required in an FMM call for those interactions that are too close to be approximated by expansions. It is also useful for debugging and testing the accuracy of the FMM call. Overloading for ::Gravitational systems, we have:

function FastMultipole.direct!(target_system, target_index, ::DerivativesSwitch{PS,VS,GS}, source_system::Gravitational, source_buffer, source_index) where {PS,VS,GS}
    @inbounds for i_source in source_index
        source_x, source_y, source_z = FastMultipole.get_position(source_buffer, i_source)
        source_strength = FastMultipole.get_strength(source_buffer, source_system, i_source)[1]
        @inbounds for j_target in target_index
            target_x, target_y, target_z = FastMultipole.get_position(target_system, j_target)
            dx = target_x - source_x
            dy = target_y - source_y
            dz = target_z - source_z
            r2 = dx*dx + dy*dy + dz*dz
            if r2 > 0
                r = sqrt(r2)
                if PS
                    dϕ = source_strength / r * FastMultipole.ONE_OVER_4π
                    FastMultipole.set_scalar_potential!(target_system, j_target, dϕ)
                end
                if VS
                    dF = SVector{3}(dx,dy,dz) * source_strength / (r2 * r) * FastMultipole.ONE_OVER_4π
                    FastMultipole.set_gradient!(target_system, j_target, dF)
                end
            end
        end
    end
end

Note that the velocity gradient is not calculated in this function. If the velocity gradient is requested, the contribution due to ::Gravitational systems will then be zero.

Running the FMM

Now that we have defined the interface functions, we can run the FMM on our Gravitational system to obtain the gravitational acceleration at each body. The fmm! function is used to perform the FMM call as:

using Random

# create system
function generate_gravitational(seed, n_bodies; radius_factor=0.1, strength_scale=1/n_bodies, bodies_fun=(x)->x)
    # random seed
    Random.seed!(seed)

    # initialize bodies
    bodies = rand(5, n_bodies)

    # scale radius
    bodies[4,:] ./= (n_bodies^(1/3)*2)
    bodies[4,:] .*= radius_factor

    # scale strength
    bodies[5,:] .*= strength_scale

    # user-defined function for arbitrary transformation
    bodies_fun(bodies)

    # create system
    return Gravitational(bodies)
end

n_bodies, rand_seed = 5_000, 123
system = generate_gravitational(rand_seed, n_bodies)

# run FMM
fmm!(system; scalar_potential=false, gradient=true, hessian=true)

println("gravitational acceleration:\n", system.potential[5:7,1:10], "...")

gravitational acceleration:
[0.02952639011390768 -0.02069619176239785 0.025308790268953336 0.05962204518034472 0.0504774301877218 -0.035587328063329385 0.07821351262466984 -0.006666310179674633 0.026898010956689013 0.05271241254316844; 0.005306546782555571 -0.017169549499724216 -0.034257888440946294 0.014788007066020278 -0.01087048265953569 0.07495477398739561 0.04177523139954798 0.0745141881656066 0.08814587758750614 0.007071761637510552; 0.026157029327728353 0.07837748716099167 -0.07404854130606213 -0.024452264548083243 -0.042227331058716004 -0.02623074799497789 -0.01386135371171003 0.03345855691484015 -0.007767314674050915 -0.02552792181490425]...

The scalar_potential and gradient arguments are set to false and true, respectively, to indicate that we want to compute the vector field (i.e. the gravitational field lines that translate to force and their gradient) but not the scalar potential.

Preallocating the Buffers

The fmm! function incurs a small overhead for allocating the buffer matrices. If our system is large and the FMM will be called more than once, we can preallocate the buffers to avoid this overhead. This is done by returning a cache of preallocated memory the first time fmm! is called, and then adding it as an additional argument to fmm! in subsequent calls as follows:

# allocate buffer cache
allocs_1 = @allocated _, cache, _ = fmm!(system; scalar_potential=false, gradient=true, hessian=true)

# run FMM with preallocated buffers
allocs_2 = @allocated fmm!(system, cache; scalar_potential=false, gradient=true, hessian=true)

println("memory allocated without the cache:  ", allocs_1, " bytes")
println("memory allocated with the cache:     ", allocs_2, " bytes")

memory allocated without the cache:  8407560 bytes
memory allocated with the cache:     7257304 bytes

Note that the second call to fmm! allocates less memory.