Multiple Systems

The most advanced features of FastMultipole include:

• the ability to solve the $n$-body problem for multiple systems simultaneously in a single FMM call
• satisfy an error tolerance by dynamically adjusting the expansion order of multipole-to-local transformations
• predict the full array of FMM tuning parameters such that an error tolerance is met

In the this section, and continuing in Automated Tuning, we will first describe how to run the FMM for multi-system problems. Then, we will describe how to impose an error tolerance. Finally, we will show how to automatically tune the FMM parameters.

Simultaneous Computation of Several Target and Source Systems

FastMutipole allows FMM to be performed efficiently on distinct target and source systems. This is done by creating two octrees: one for sources, and one for targets [3]. For example, say we desire to know the influence of one system of point masses on another, but not vice versa. This is done by passing both systems to the fmm! function with the target system first, as:

using FastMultipole
using Random # needed for `gravitational.jl`

gravitational_path = normpath(joinpath(splitdir(pathof(FastMultipole))[1], "..", "test", "gravitational.jl"))
include(gravitational_path)

target_system = generate_gravitational(123, 1000)
source_system = generate_gravitational(321, 1000)

fmm!(target_system, source_system)
println("gravitational acceleration: ", target_system.potential[5:7,1:10], "...")

gravitational acceleration: [0.08493884822133471 -0.11084372433205175 0.041347632985565604 -0.1863621517367257 -0.08011810805054259 0.04306673245903067 0.0026037448998806746 0.03587428613837394 -0.08002640294442173 0.15009451094982013; 0.002073840920276872 0.3215930836528468 -0.03967273201136324 -0.005960953149224658 0.03828114050071747 0.14489807182906364 0.010013929319383668 0.05641013581893285 0.021727973779433876 0.007681084388235132; 0.023957945365164404 -0.09630611059971712 -0.025023566343257052 0.1337146690539532 -0.04538444503621948 -0.07362887885237299 0.08060484661663683 -0.015529897358683238 -0.015618344404806379 0.10559335815763894]...

In practice, the source system might be a collection of systems, composed of a variety of datastructures. So long as the interface functions are defined for each system, we can pass a tuple of any number of systems as the source and or target. For example, we could evaluate the influence of a system of point masses, point vortices, and vortex filaments on two different target systems:

using LinearAlgebra

# include vortex filament and particle models and interface functions
vortex_path = normpath(joinpath(splitdir(pathof(FastMultipole))[1], "..", "test", "vortex.jl"))
include(vortex_path)
filament_path = normpath(joinpath(splitdir(pathof(FastMultipole))[1], "..", "test", "vortex_filament.jl"))
include(filament_path)

# generate systems
n_bodies = 2000
target_one = generate_gravitational(123, n_bodies)
target_two = generate_vortex(124, n_bodies)
source_one = generate_gravitational(125, n_bodies)
source_two = generate_vortex(126, n_bodies)
source_three = generate_filament_field(n_bodies, n_bodies^0.333, 127; strength_scale=1/n_bodies)

# run FMM
fmm!((target_one, target_two), (source_one, source_two, source_three);
    scalar_potential=false, gradient=true, hessian=(false, true))

v1 = target_one.potential[5:7,:]
v2 = target_two.gradient_stretching[1:3,:]

# run direct
target_one.potential .= 0.0
target_two.potential .= 0.0
target_two.gradient_stretching .= 0.0
direct!((target_one, target_two), (source_one, source_two, source_three); scalar_potential=false, gradient=true, hessian=(false, true))

# test accuracy
println("max error in target one: ", maximum(abs.(target_one.potential[5:7,:] .- v1)))
println("max error in target two: ", maximum(abs.(target_two.gradient_stretching[1:3,:] .- v2)))

max error in target one: 1.4007503612439853e-5
max error in target two: 1.0289615966441179e-5

Note that scalar_potential, gradient, and hessian can be passed as a single boolean or as a tuple of booleans, one for each target system. This allows the user to specify which values are desired for each target system, and avoids unnecessary calculations for values that are not needed. In this case, we have set scalar_potential=false and gradient=true to indicate all target systems, but a tuple hessian=(false,true) to indicate different settings for each. It is worth remembering that these switches must be implemented by the user when overloading the direct! function for each system to act as a source.