Predefined kernels

While it is certainly possible to define your own kernel, doing so efficiently can be non-trivial and time-consuming. IFGF.jl provides a few predefined kernels that are commonly used in practice. If you need a kernel that is not listed here, or if you have implemented a custom kernel and would like for it to be added to the library, please open an issue and/or PR on the GitHub repository.

The kernel-specific API is heavily inspired by the FMM3D library, and provides functions to efficiently compute sums of the form

\[u(\boldsymbol{x}_i) = \sum_{j=1}^N G(\boldsymbol{x}_i - \boldsymbol{y}_j) c_j +\gamma_{1,\boldsymbol{y}} G(\boldsymbol{x}_i - \boldsymbol{y}_j) v_j,\]

where $\boldsymbol{x_i}$ are the target locations, $\boldsymbol{y}_j$ are the source locations, $G(\boldsymbol{x} - \boldsymbol{y})$ is (up to a constant factor) the fundamental solution of a given partial differential operator, $\gamma_{1,\boldsymbol{y}}$ its (generalized) Neumann trace, and $c_i,v_i$ are input vectors of the appropriate size. IFGF.jl computes the sum above in log-linear complexity.

Each predefined kernel comes with an associated plan function that precomputes the necessary data structure required to efficiently apply the forward map. If you need to compute the forward map multiple times with the same kernel, it is a good idea to build a plan.

In the following, we provide a brief overview of the predefined kernels.

Single-precision

All predefined kernels support single-precision arithmetic. If you do not need double-precision, you can save memory and computation time by using single-precision to represent your data.

Laplace 3D operator

IFGF.laplace3d — Function

laplace3d(targets, sources; charges = nothing, dipvecs = nothing, grad =
false, kwargs...)

Compute the sum

\[u(\boldsymbol{x}_i) = \sum_{j=1}^n \frac{c_j}{|| \boldsymbol{x}_i - \boldsymbol{y}_j ||} - \frac{(\boldsymbol{x}_i - \boldsymbol{y}_j)^\top v_j}{|| \boldsymbol{x}_i - \boldsymbol{y}_j ||^3} ,\]

where $c_j$ are the charges and v_j are the dipvecs.

Input:

targets::Matrix{T}: 3 x n matrix of target points.
sources::Matrix{T}: 3 x m matrix of source points.
charges::Vector{T}: vector of n charges.
dipvecs::Matrix{T}: 3 x m matrix of dipole vectors.
grad::Bool: if true, the gradient is computed instead
kwargs...: additional keyword arguments passed to assemble_ifgf (e.g. tol)

Output

u: the potential at the target points if grad=false, or the gradient as a 3 x n matrix if grad=true.

source

Forward map:

using IFGF
m = 100_000
targets = sources = rand(3, m)
charges = rand(m)
dipvecs = rand(3,m)
IFGF.laplace3d(targets, sources; charges, tol=1e-4)

100000-element Vector{Float64}:
  97275.90663479632
  85615.07647178617
 108228.48388642876
  82272.8660121476
  68054.01977632476
  99911.24126839136
  91588.6463015972
  84028.6886121075
  90198.04368304361
 102134.2694682958
      ⋮
  79445.94879050551
 101466.92017509496
  91116.8007449042
  91815.65869005359
  80944.18713845665
  78168.26765462913
  78948.10715639617
  83605.71958075691
  88114.51117406323

Plan and forward map:

L = IFGF.plan_laplace3d(targets, sources; charges, dipvecs, tol=1e-4)
out  = zero(charges)
IFGF.laplace3d!(out, L; charges, dipvecs)

100000-element Vector{Float64}:
  -18793.92360803263
   43952.618943565394
  207207.83046361338
 -100803.1205071884
  114864.82713589253
  138073.27412826527
   81247.419487003
  131929.20458311465
   96216.5459758685
  187471.019137077
       ⋮
  125941.57101007292
   88955.8471379106
  177587.80608160037
  -59029.41660573503
   38682.625442183125
   95779.01120395274
 -104871.95295999388
   89380.29061228389
  140543.6924905796

Helmholtz 3D operator

IFGF.helmholtz3d — Function

helmholtz3d(k, targets, sources; charges = nothing, dipvecs = nothing, grad =
false, kwargs...)

Compute the sum

\[u(\boldsymbol{x}_i) = \sum_{j=1}^n \frac{c_j e^{i k || \boldsymbol{x}_i - \boldsymbol{y}_j ||}}{|| \boldsymbol{x}_i - \boldsymbol{y}_j ||} + \nabla_{\boldsymbol{y}}^\top \left( \frac{e^{i k || \boldsymbol{x}_i - \boldsymbol{y}_j ||}}{|| \boldsymbol{x}_i - \boldsymbol{y}_j ||} \right) v_j,\]

where $c_j$ are the charges and v_j are the dipvecs.

Input:

k: the wavenumber.
targets::Matrix{T}: 3 x n matrix of target points.
sources::Matrix{T}: 3 x m matrix of source points.
charges::Vector{Complex{T}}: vector of n charges.
dipvecs::Matrix{Complex{T}}: 3 x m matrix of dipole vectors.
grad::Bool: if true, the gradient is computed instead
kwargs...: additional keyword arguments passed to assemble_ifgf.

Output

u: the potential at the target points if grad=false, or the gradient as a 3 x n matrix if grad=true.

source

Forward map:

using IFGF
m = 100_000
targets = sources = rand(3, m)
charges = rand(ComplexF64,m)
dipvecs = rand(ComplexF64,3,m)
k = 2π
IFGF.helmholtz3d(k, targets, sources; charges, tol=1e-4)

100000-element Vector{ComplexF64}:
  -3243.300512356839 - 8117.369304174505im
  -56933.31907836979 - 24780.380794319313im
 -62433.528917150965 - 26728.817007239188im
    4917.94479504552 - 5143.59908848029im
 -11009.092685643156 - 15854.40827673315im
  -37709.83121433167 - 19906.23743634126im
  -35961.74275210078 - 19805.578260938953im
 -25877.739606858253 - 17717.391802810926im
  -48518.00626509611 - 24589.413318434665im
 -22456.313328955544 - 13955.025365157679im
                     ⋮
 -35071.345463282036 - 20703.336984143127im
  -5366.507578187988 - 7548.747491729477im
  -4675.841885738251 - 12612.18123993253im
 -4011.2166413238892 - 10426.576890068085im
 -51711.063775224466 - 22225.17066957883im
 -32116.903793523947 - 18578.3913347249im
 -16482.333140847226 - 14279.837444868946im
   4498.781930463909 - 4057.2636284080145im
 -31856.046873321488 - 17828.690694694244im

Plan and forward map:

L = IFGF.plan_helmholtz3d(k, targets, sources; charges, dipvecs, tol=1e-6)
out  = zero(charges)
IFGF.helmholtz3d!(out, L; charges, dipvecs)

100000-element Vector{ComplexF64}:
  -106769.2431649079 - 27799.981683397986im
  -409542.6130361583 - 100846.55408547701im
  248748.67226259015 + 46559.7464700741im
 -22871.373174728975 - 45955.35642446423im
 -232002.27382100132 - 92748.59900914202im
  241443.42639861532 + 54856.870487335706im
  234477.19374100555 + 46509.008244638244im
  3705.6447638072204 + 28516.18538331025im
 -291769.36395134847 - 56535.53549137121im
 -126878.12655769344 - 44389.99305956374im
                     ⋮
   112098.1876399829 - 4246.877042546646im
 -265676.78068018804 - 125761.84021822008im
  174927.91877509234 + 64912.58601866248im
  -39410.71451174546 - 59822.752792897736im
 -102871.57495747026 - 43938.435761248315im
   95373.32737615955 + 11222.300393959598im
 -18854.007963548713 + 12469.881253577892im
 -44142.187357732764 - 110295.63230902613im
  -399819.0881164571 - 131495.57635338441im

Stokes 3D operator

Warning

Stokes 3D operator should be considered experimental. Please report any issues you encounter.

IFGF.stokes3d — Function

stokes3d(targets, sources; stoklet = nothing, strslet = nothing, strsvec =
nothing, grad = false, kwargs...)

Compute the sum

\[u(\boldsymbol{x}_i) = \sum_{j=1}^n \left( \frac{I}{d} + \frac{\boldsymbol{r}_{ij}\boldsymbol{r}_{ij}^\top}{2d_{ij}^3} \right) \boldsymbol{\sigma}_j + \frac{3(\boldsymbol{\nu}_j \cdot \boldsymbol{r}_{ij})(\boldsymbol{\mu}_j\cdot\boldsymbol{r}_{ij}) \boldsymbol{r}_{ij}}{d_{ij}^5},\]

where $\boldsymbol{r}_{ij} = |\boldsymbol{x}_i - \boldsymbol{y}_j|$, $d_{ij} = |\boldsymbol{r_{ij}}|$, $\boldsymbol{\sigma}_j$ are the Stokeslet strengths, $\boldsymbol{\mu}_j$ are the stresslet strengths, and $\boldsymbol{\nu}_j$ are the stresslet orientations.

Input:

targets::Matrix{T}: 3 x n matrix of target points.
sources::Matrix{T}: 3 x m matrix of source points.
stoklet::Matrix{T}: 3 × n matrix of stokeslets
strslet::Matrix{T}: 6 x m matrix of stresslet densities (indices 1:3) and orientations (indices 4:6)
grad::Bool: if true, the gradient is computed instead
kwargs...: additional keyword arguments passed to assemble_ifgf.

Output

u: 3 x n matrix giving the velocity at the target points if grad=false, or 3 x 3 × n matrix representing the gradient of u if grad=true.

source

Forward map:

using IFGF
m = 100_000
targets = sources = rand(3, m)
stoklet  = rand(3,m)
strslet = nothing
IFGF.stokes3d(targets, sources; stoklet, strslet, tol=1e-4)

3×100000 Matrix{Float64}:
 76455.2  55198.2  51014.8  58654.1  …  66863.1  56297.1  54452.2  74033.3
 76444.0  42446.6  60076.9  62257.6     58679.4  55132.3  63695.4  75472.6
 74170.0  52306.3  59855.9  63533.5     68993.6  56325.7  63275.4  75223.5

Plan and forward map:

L = IFGF.plan_stokes3d(targets, sources; stoklet, strslet, tol=1e-4)
out = zero(stoklet)
IFGF.stokes3d!(out, L; stoklet)

3×100000 Matrix{Float64}:
 76455.2  55198.2  51014.8  58654.1  …  66863.1  56297.1  54452.2  74033.3
 76444.0  42446.6  60076.9  62257.6     58679.4  55132.3  63695.4  75472.6
 74170.0  52306.3  59855.9  63533.5     68993.6  56325.7  63275.4  75223.5