BiotSavart

BiotSavart

Module for estimation of Biot–Savart integrals along vortex filaments using fast Ewald splitting.

ParamsBiotSavart{T <: AbstractFloat}

Contains parameters for calculation of Biot–Savart integrals using fast Ewald splitting.

The type parameter T corresponds to the precision used in computations (typically Float64 or Float32).

Construction

ParamsBiotSavart([T = Float64]; Γ, a, α, Ls, Ns, rcut, optional_kws...)

where the optional parameter T sets the numerical precision.

Mandatory and optional keyword arguments are detailed in the extended help below.

See also BiotSavart.autotune for an alternative way of setting Biot–Savart parameters.

Extended help

Mandatory keyword arguments

• Γ::Real: vortex circulation (assumed constant);

• a::Real: vortex core size (assumed constant);

• α::Real: Ewald splitting parameter (inverse length scale). One can set α = Zero() to efficiently disable long-range computations.

• Ls::Union{Real, NTuple{3, Real}}: domain period in each Cartesian direction. If a single value is passed (e.g. Ls = 2π), it is assumed that periods are the same in all directions.

  One can set Ls = Infinity() to disable periodicity. This should be done in combination with α = Zero().

• Ns::Dims{3}: dimensions of physical grid used for long-range interactions. This parameter is not required if α = Zero().

• rcut: cutoff distance for computation of short-range interactions. For performance and practical reasons, the cutoff distance must be less than half the cell unit size in each direction, i.e. rcut < minimum(Ls) / 2. This parameter is not required if α = Zero().

Optional keyword arguments (and their defaults)

General

• quadrature::StaticSizeQuadrature = GaussLegendre(3): quadrature rule for short- and long-range interactions. For example, if quadrature = GaussLegendre(4), then 4 evaluations of the Biot–Savart integrand will be done for each filament segment.

Short-range interactions

• backend_short::ShortRangeBackend: backend used to compute short-range interactions. The default is CellListsBackend(2), unless periodicity is disabled, in which case NaiveShortRangeBackend() is used. See ShortRangeBackend for a list of possible backends;

Long-range interactions

• backend_long::LongRangeBackend = NonuniformFFTsBackend(): backend used to compute long-range interactions. See LongRangeBackend for a list of possible backends;

• longrange_truncate_spherical = false: if true, perform a spherical truncation in Fourier space, discarding all wavenumbers such that $|\bm{k}| > k_{\text{max}}$. This is not recommended as it leads to precision loss, and should be used for testing only (in particular, for verifying error estimates which assume this kind of truncation).

Local self-induced velocity

• Δ = 0.25: coefficient appearing in the local self-induced velocity (LIA term), which depends on the vorticity profile at the vortex core.

  Some common values of Δ are:

  • Δ = 0.25 for a constant vorticity profile (default);

  • Δ = 0.5 for a hollow vortex;

  • Δ ≈ 0.905 ≈ 0.558 + ln(2) / 2 for a Gaussian vorticity profile (taking a as the Gaussian standard deviation σ);

  • Δ ≈ 0.615 for a Gross–Pitaevskii vortex with healing length a.

  See Saffman (1992), sections 10.2–10.3 for the first three.

• lia_segment_fraction = nothing: can be used to indicate that the LIA term should be evaluated over a fraction of the two segments surrounding a node. In this case, it should be a real value in $(0, 1]$. The default (nothing) is equivalent to 1, and means that the LIA term is evaluated over the full segments. If smaller than 1, the velocity induced by the excluded part of the segments will be evaluated using the regular Biot–Savart law (using quadratures within each subsegment). This may improve accuracy, especially when the discretisation distance is relatively large. Since this means integrating near the singularity of the BS integral, this integral is by default estimated using adaptive quadratures (see quadrature_near_singularity below).

• quadrature_near_singularity = AdaptiveTanhSinh(T; nlevels = 5): quadrature rule to be used when integrating near a singularity, in particular when lia_segment_fraction is enabled. By default an adaptive quadrature rule AdaptiveTanhSinh is used, which is generally accurate but can be costly. One can also pass a StaticSizeQuadrature such as GaussLegendre, but in that case accuracy is not guaranteed.

BiotSavart.autotune(fs::AbstractVector{<:AbstractFilament{T}}, [β::Real = 3.5]; kwargs...) -> ParamsBiotSavart{T}

Generate Biot-Savart parameters optimised for a given filament geometry and accuracy parameter $β$.

Based on the given value of $β$, this function will try to automatically set the following Ewald-related parameters:

• inverse splitting parameter $α$;
• cut-off distance $r_{\text{cut}} = β / α$;
• Fourier-space cut-off $k_{\text{max}} = 2αβ ≈ πM/L$, where $M$ is the grid size and $L$ the domain period.

In practice, this function will try to find the value of $α$ which minimises the computation of the velocity of the filaments fs.

See also ParamsBiotSavart.

Extended help

Mandatory keyword arguments

• Γ::Real: vortex circulation (assumed constant);

• a::Real: vortex core size (assumed constant);

• Ls::Union{Real, NTuple{3, Real}}: domain period in each Cartesian direction. If a single value is passed (e.g. Ls = 2π), it is assumed that periods are the same in all directions.

Optional keyword arguments

This function accepts the same keyword arguments as ParamsBiotSavart and with the same default values. In particular, the default short- and long-range backends are:

• backend_short::ShortRangeBackend = CellListsBackend(2);

• backend_long::LongRangeBackend = NonuniformFFTsBackend(σ = 1.5, m = HalfSupport(4)).

Autotuning parameters

The following keyword arguments can be used to control autotuning:

• nruns = 4: number of Biot–Savart calculations per value of $α$. The minimum elapsed time among all runs will be used in autotuning;

• Cstart: initial guess for non-dimensional factor $C$ (see Autotuning algorithm below). Default values are 1.5 (pure CPU) and 4.0 (CPU + GPU).

• ΔC = 0.1: increment of non-dimensional factor $C$;

• verbose = false: if true, print autotuning information.

Autotuning algorithm

The autotuning algorithm basically consists in trying different values of $α$, which we write under the form:

\[α = C \, {( N / V )}^{1/3}\]

where $N$ is the total number of filament nodes and $V$ the domain volume. The parameter that is varied is the non-dimensional factor $C$.

For now the algorithm is quite basic. We try different values around C = Cstart using increments of ΔC. The parameters giving the fastest runtime are returned.

Typical values of $β$ and NUFFT parameters

The following table roughly relates accuracy (in number of digits) and values of $β$, as detailed in Polanco (2024):

The last column is the size of the NUFFT half-support $w$ which ensures sufficient NUFFT accuracy. The given values assume a NUFFT oversampling factor $σ = 1.5$ and a (backwards) Kaiser–Bessel spreading kernel, which are the default when using the NonuniformFFTsBackend. Currently, this parameter is not automatically selected by this function. In other words, knowing the required value of $w$, one can pass:

backend_long = NonuniformFFTsBackend(m = HalfSupport(w))

as a keyword argument to this function.

BiotSavart.circulation(p::ParamsBiotSavart) -> Γ

Return the circulation Γ associated to each vortex.

BiotSavart.periods(p::ParamsBiotSavart) -> (Lx, Ly, Lz)

Return the domain periods in each direction.

BiotSavart.domain_is_periodic(p::ParamsBiotSavart) -> Bool
BiotSavart.domain_is_periodic(Ls::NTuple) -> Bool

Check whether the domain is periodic.

Returns true if the domain is periodic in all directions, false otherwise.

init_cache(
    p::ParamsBiotSavart, fs::AbstractVector{<:AbstractFilament};
    timer = TimerOutput("BiotSavart"),
) -> BiotSavartCache

Initialise caches for computing Biot–Savart integrals.

velocity_on_nodes!(
    vs::AbstractVector{<:AbstractVector{<:Vec3}},
    cache::BiotSavartCache,
    fs::AbstractVector{<:AbstractFilament};
    kws...,
) -> vs

Compute velocity induced by vortex filaments on filament nodes.

Velocities induced by vortex filaments fs are written to vs.

This is the same as calling compute_on_nodes! when only the velocity is needed.

Usually, fs is a vector containing all the vortex filaments in the system. In that case, vs must be a vector of vectors, which will contain the velocities of all filament nodes. The length of vs[i] must be equal to the number of nodes in the filament fs[i].

The vector of velocities where the output will be written may be initialised using one of the following lines (all are exactly equivalent):

vs = map(similar ∘ nodes, fs)
vs = [similar(nodes(f)) for f ∈ fs]
vs = similar.(nodes.(fs))

which initialise a velocity vector for each node of each filament (see also nodes).

See compute_on_nodes! for a list of accepted keyword arguments.

compute_on_nodes!(
    fields::NamedTuple{Names, NTuple{N, V}},
    cache::BiotSavartCache,
    fs::AbstractVector{<:AbstractFilament};
    LIA = Val(true),
    shortrange = true,
    longrange = true,
    callback_vorticity = identity,
) where {Names, N, V <: AbstractVector{<:VectorOfVec}}

Compute velocity and/or streamfunction on filament nodes.

The first argument contains one or more output fields to compute. It is usually of length 1 or 2, and can contain fields named velocity and streamfunction.

For example, to compute both velocity and streamfunction on the nodes of filaments fs:

# Initialise fields to compute (vectors of vectors)
vs = map(similar ∘ nodes, fs)  # one velocity vector per filament node
ψs = map(similar, vs)

# The first argument to `compute_on_nodes!` must have the following form.
# One can also choose to pass just one of the two fields.
fields = (;
    velocity = vs,
    streamfunction = ψs,
)

cache = BiotSavart.init_cache(...)
compute_on_nodes!(fields, cache, fs)

Extended help

Disabling local terms / computing only local terms

One may disable computation of the locally-induced velocity and streamfunction (LIA term) by passing LIA = Val(false). Conversely, one can pass LIA = Val(:only) to compute only the LIA term. This can be useful for splitting the induced filament velocities/streamfunctions onto local and non-local parts.

Disabling short-range or long-range interactions

It is also possible to disable the short-range or long-range component of Ewald splitting, if one only wants to compute one of the two components. To do this, pass either shortrange = false or longrange = false.

Note that the short-range component includes the local (LIA) term as well as the short-range correction term for long-range interactions. Therefore, setting shortrange = false disables both these terms.

Accessing vorticity in Fourier space

The computation of long-range quantities involves estimating the Fourier coefficients of the vorticity field associated to the vortex filaments. These coefficients are truncated to some maximum wavenumber $k_{\text{max}}$ in each Cartesian direction. This information can be useful for other things, for instance computing energy spectra.

One can use the callback_vorticity argument to access the vorticity in Fourier space, before it is replaced by the coefficients of streamfunction and/or velocity. This argument should be a function callback_vorticity(cache) which takes a LongRangeCache. The callback should not modify anything inside the cache, or otherwise the streamfunction and velocity computed by this function will likely be wrong. Moreover, one can call get_longrange_field_fourier to get the vorticity field from within the callback. Of course, this callback never be called if long-range computations are disabled.

Note that, when the callback is called, the vorticity coefficients returned by get_longrange_field_fourier don't have the right physical dimensions as they have not yet been multiplied by $Γ/V$ (where $V$ is the volume of a unit cell). Note that $Γ/V$ is directly available in cache.common.ewald_prefactor. Besides, the vorticity coefficients at this stage have not yet been Gaussian-filtered according to Ewald's method.

An example of how to compute the (large-scale) kinetic energy associated to the Fourier-truncated vorticity field:

using Adapt: adapt  # useful in case FFTs are computed on the GPU

E_from_vorticity = Ref(0.0)  # "global" variable updated when calling compute_on_nodes!

function callback_vorticity(cache::LongRangeCache)
    (; ewald_prefactor,) = cache.common
    (; field, wavenumbers, state,) = BiotSavart.get_longrange_field_fourier(cache)
    @assert state.quantity == :vorticity
    @assert state.smoothed == false
    # To make things simple, we copy data to the CPU if it's on the GPU.
    wavenumbers = adapt(Array, wavenumbers)
    uhat = adapt(Array, field)::NTuple{3}  # (ωx, ωy, ωz) in Fourier space
    with_hermitian_symmetry = BiotSavart.has_real_to_complex(cache)  # this depends on the long-range backend
    @assert with_hermitian_symmetry == wavenumbers[1][end] > 0
    γ² = ewald_prefactor^2  # = (Γ/V)^2 [prefactor not included in the vorticity]
    E = 0.0
    for I ∈ CartesianIndices(uhat[1])
        k⃗ = map(getindex, wavenumbers, Tuple(I))
        kx = k⃗[1]
        factor = (!with_hermitian_symmetry || kx == 0) ? 0.5 : 1.0
        k² = sum(abs2, k⃗)
        if !iszero(k²)
            ω⃗ = (uhat[1][I], uhat[2][I], uhat[3][I])  # Fourier coefficient of the vorticity
            E += γ² * factor * sum(abs2, ω⃗) / k²
        end
    end
    E_from_vorticity[] = E  # update value of "global" variable
    nothing
end

compute_on_nodes!(...; callback_vorticity)

BiotSavart.kelvin_wave_period(p::ParamsBiotSavart, λ::Real) -> Real

Return the period $T(λ)$ associated to Kelvin waves of wavelength $λ$.

This can be convenient for setting the timestep dt associated to a filament discretisation distance δ. The timestep should typically be proportional to the period of the Kelvin waves of wavelength δ.

The Kelvin wave period is $T(λ) = 2π/ω(k)$ where $k = 2π/λ$ is the wavenumber associated to $λ$ and $ω(k)$ is the Kelvin wave dispersion relation:

\[ω(k) = \frac{Γ k^2}{4π} \left[ \ln\left( \frac{2}{k a} \right) - γ + \frac{1}{2} - Δ \right]\]

where $γ ≈ 0.5772$ is the Euler–Mascheroni constant.

ShortRangeBackend

Abstract type denoting the backend used for computing short-range interactions.

Implemented backends

• CellListsBackend: most efficient when the cutoff radius is much smaller than the domain size. Can only be used with periodic boundary conditions.

• NaiveShortRangeBackend: usually less efficient as it needs to compute distances between all pairs of locations.

Extended help

Implementation details

A BACKEND <: ShortRangeBackend must implement the function:

init_cache_short(c::ParamsCommon, p::ParamsShortRange{<:BACKEND}, fs::AbstractVector{<:AbstractFilament}, to::TimerOutput)

which should return a ShortRangeCache.

It may also implement the function max_cutoff_distance.

NaiveShortRangeBackend <: ShortRangeBackend

Naive computation of short-range interactions.

CellListsBackend <: ShortRangeBackend
CellListsBackend(nsubdiv::Int = 1)

Compute short-range interactions using the cell lists algorithm.

This backend can be significantly faster than the NaiveShortRangeBackend when the cutoff distance r_cut is much smaller than the domain period L (roughly when r_cut ≲ L / 10).

Optionally, one can choose to subdivide each cell (of size ≈ r_cut) onto nsubdiv subcells. In practice, a value of 2 or 3 can significantly improve performance compared to no subdivision (1).

Note that, with this backend, the cutoff distance must satisfy r_cut ≤ M / (2M + 1) * L where M = nsubdiv.

This backend does not support non-periodic domains.

See PeriodicCellList and Wikipedia for more details.

ShortRangeCache

Abstract type describing the storage of data required to compute short-range interactions.

The init_cache_short function returns a concrete instance of a ShortRangeCache.

Interface

Fields

The following fields must be included in a cache:

• params :: ParamsShortRange parameters for short-range computations;

• to :: TimerOutput for measuring time spent on different functions.

max_cutoff_distance(::ShortRangeBackend, L::Real) -> r
max_cutoff_distance(::ShortRangeBackend, Ls::NTuple{3, Real}) -> r

Return the maximum cut-off distance r_cut allowed by the backend for a given domain period L.

This is usually close to L/2, but the actual value can depend on implementation details of the backend. For example, the CellListsBackend requires a slightly smaller distance, in the range L/3 ≤ r_max < L/2 depending on the backend parameters.

init_cache_short(
    pc::ParamsCommon, p::ParamsShortRange,
    fs::AbstractVector{<:AbstractFilament},
    to::TimerOutput,
) -> ShortRangeCache

Initialise the cache for the short-range backend defined in p.

process_point_charges!(c::ShortRangeCache, data::PointData)

Process list of point charges.

This is useful for short-range backends like CellListsBackend, which needs to

Must be called after add_point_charges! and before computing any short-range quantities (using add_short_range_fields!).

add_short_range_fields!(
    fields::NamedTuple{Names, NTuple{N, V}},
    cache::ShortRangeCache,
    f::AbstractFilament;
    LIA = Val(true),
)

Compute short-range Biot-Savart integrals.

Adds the results onto fields. See compute_on_nodes! for more details.

Setting LIA = Val(false) allows to disable computation of the localised induction approximation (LIA) term. In that case, that term should be computed separately using local_self_induced.

Before calling this function, one must first call add_point_charges! and then process_point_charges!.

local_self_induced_velocity(
    f::AbstractFilament, i::Int, [prefactor::Real];
    a::Real, [Γ::Real],
    Δ = 0.25, quad = nothing, fit_circle = false,
    segment_fraction = nothing,
)

Compute local self-induced velocity of filament node f[i].

This corresponds to the LIA term (localised induction approximation).

This is the same as local_self_induced(Velocity(), f, i, [prefactor]; kwargs...). See also local_self_induced.

local_self_induced(
    q::OutputField, f::AbstractFilament, i::Int, [prefactor::Real];
    a::Real, Δ::Real = 0.25, segment_fraction = nothing, quad = nothing, [Γ],
)

Compute localised induction approximation (LIA) term at node f[i].

Possible output fields to be passed as first argument are Velocity() and Streamfunction().

One must pass either the circulation Γ as a keyword argument, or a precomputed prefactor which is usually equal to Γ / 4π (both for velocity and streamfunction).

Mandatory arguments

• the vortex core size a as a keyword argument;

• either the vortex circulation Γ as a keyword argument, or the precomputed prefactor Γ / 4π as a positional argument.

Optional arguments

• the optional parameter Δ sets the effect of the vorticity profile (see ParamsBiotSavart for details);

• the optional parameter segment_fraction can be used to indicate that the LIA term should only be computed over a fraction of the segments neighbouring the node f[i]. In that case, it should be a real number in $(0, 1]$ (where 1 corresponds to integrating over the full segments, which the default);

• a quadrature rule can be passed via quad, which can improve the estimation of the LIA term and the stability of the solver (even a 1-point quadrature rule can importantly improve stability!);

• if quad = nothing, one may set fit_circle = true to estimate the binormal vector by fitting a circle passing through 3 neighbouring nodes (as done in Schwarz PRB 1985), instead of using local derivatives. For now, this is only implemented for q = Velocity(), and it may be removed in the future.

nearby_charges(c::ShortRangeCache, x⃗::Vec3)

Return an iterator over the charges that are "close" to the location x⃗.

LongRangeBackend

Abstract type denoting the backend to use for computing long-range interactions.

Implemented backends

• NonuniformFFTsBackend: estimates long-range interactions via the non-uniform fast Fourier transform (NUFFT) using the NonuniformFFTs.jl package;

• FINUFFTBackend: estimates long-range interactions via the NUFFT using the FINUFFT library;

• CuFINUFFTBackend: estimates long-range interactions via the NUFFT using the CUDA implementation of the FINUFFT library. CUDA.jl must be loaded before using this backend (CUDA devices only);

• ExactSumBackend: computes long-range interactions using exact Fourier sums. This is really inefficient and should only be used for testing.

Extended help

Implementation details

The following functions must be implemented by a BACKEND <: LongRangeBackend:

• init_cache_long_ewald(c::ParamsCommon, p::ParamsLongRange{<:BACKEND}, to::TimerOutput) -> LongRangeCache.

• has_real_to_complex,

• expected_period (optional),

• folding_limits (optional),

• KernelAbstractions.get_backend (required for GPU-based backends).

NonuniformFFTsBackend <: LongRangeBackend

Compute long-range interactions using the NonuniformFFTs.jl package.

This backend may be faster than other NUFFT-based backends since it allows real valued non-uniform data, meaning that we can use real-to-complex FFTs to accelerate computations.

Transforms can be performed either on the CPU (parallelised with threads, default) or on a single GPU (in principle any kind of GPU should work, but only CUDA has been tested). This must be set via the device argument (see below).

Optional arguments

The signature of NonuniformFFTsBackend is:

NonuniformFFTsBackend([device = CPU()]; σ = 1.5, m = HalfSupport(4), kws...)

where all arguments are passed to NonuniformFFTs.jl.

Using a GPU

Transforms are run on all available CPUs by default. To use a GPU, pass the corresponding KernelAbstractions.jl backend as the only positional argument. For example, to use a CUDA device:

using CUDA
backend_long = NonuniformFFTsBackend(CUDABackend(); kwargs...)

On AMD GPUs the following should work:

using AMDGPU
backend_long = NonuniformFFTsBackend(ROCBackend(); kwargs...)

Keyword arguments

Some relevant keyword arguments are:

• σ = 1.5: upsampling factor, which must be larger than 1. Usual values are between 1.25 (smaller FFTs, less accurate) and 2.0 (larger FFTs, more accurate). Other values such as 1.5 (default) also work;

• m = HalfSupport(4): the half-width of the NUFFT kernels. Larger means higher accuracy;

• fftw_flags = FFTW.MEASURE: flags passed to the FFTW planner (ignored on GPU devices).

The default parameters (σ = 1.5, m = HalfSupport(4)) correspond to a relative NUFFT tolerance of $∼10^{-6}$.

See the NonuniformFFTs.jl docs for a full list of possible keyword arguments.

Effect of parameters on accuracy

The following table roughly relates accuracy (in number of digits) and NUFFT parameters, as detailed in Polanco (2024):

The last column is the associated value of the accuracy parameter $β$ in Ewald's method as formulated in Polanco (2024). Once one has set $β$ and Ewald's splitting parameter $α$ (an inverse lengthscale), the cut-offs in physical and Fourier space are $r_{\text{cut}} = β / α$ and $k_{\text{max}} = 2βα$. In this formulation, $β$ controls the method accuracy while $α$ is tuned to maximise performance.

FINUFFTBackend <: LongRangeBackend

Compute long-range interactions using the FINUFFT.jl package.

To use this backend, one first needs to load FINUFFT by doing

using FINUFFT

This package provides a Julia interface to the FINUFFT C++ library, which enables efficient and accurate computation of non-uniform fast Fourier transforms (NUFFTs) based on an "exponential of a semi-circle" kernel.

Computations can be parallelised using threads.

Optional arguments

The signature of FINUFFTBackend is:

FINUFFTBackend(; tol = 1.0e-6, kws...)

where all arguments are passed to FINUFFT.

Some relevant options are:

• tol = 1.0e-6 tolerance in NUFFT computations;

• upsampfac = 1.25 upsampling factor. Must be either 1.25 (usually faster) or 2.0 (required to exceed 9 digits of accuracy);

• nthreads = Threads.nthreads() number of threads to use. By default, all threads available to Julia are used;

• fftw = FFTW.MEASURE flags passed to the FFTW planner.

Other options described in the FINUFFT docs and not listed above are also accepted.

CuFINUFFTBackend <: LongRangeBackend

GPU version of FINUFFTBackend.

To use this backend, one first needs to load both FINUFFT and CUDA by doing

using FINUFFT
using CUDA

Works with Nvidia GPUs only.

Optional arguments

The signature of CuFINUFFTBackend is:

CuFINUFFTBackend(; tol = 1.0e-6, kws...)

where all arguments are passed to cuFINUFFT.

Some relevant options are:

• tol = 1.0e-6 tolerance in NUFFT computations;

• upsampfac = 1.25 upsampling factor. Should be either 1.25 or 2.0;

• gpu_device::CUDA.CuDevice: useful if multiple GPUs are available. By default the currently active CUDA device is used, i.e. gpu_device = CUDA.device().

See the cuFINUFFT docs for details and other possible options.

ExactSumBackend <: LongRangeBackend

Compute long-range interactions "exactly" (up to filament discretisation errors) using sums of Fourier series across non-uniform points.

This should only be used for testing, as it is very slow and scales very badly with the number of non-uniform points.

BiotSavart.get_longrange_field_fourier(cache) -> NamedTuple

Obtain long-range field used to compute long-range interactions.

The input cache can be a BiotSavartCache or a LongRangeCache.

This function returns a NamedTuple with the fields:

• field: a Tuple (ux, uy, uz) describing a vector field in Fourier space. Each component is a 3D matrix of complex values.

• wavenumbers: a Tuple (kx, ky, kz) with the wavenumbers associated to the grid in Fourier space.

• state: allows to know what the returned field actually represents (vorticity, velocity, ...). See LongRangeCacheState for details.

LongRangeCache

Abstract type describing the storage of data required to compute long-range interactions.

The init_cache_long function returns a concrete instance of a LongRangeCache (or NullLongRangeCache(), if long-range computations were disabled by setting α = Zero()).

Useful fields

Most useful fields of a cache::LongRangeCache are in the cache.common field. In particular, cache.common contains the fields:

• wavenumbers_d::NTuple{3, AbstractVector}: Fourier wavenumbers in each direction;
• uhat_d::StructArray{Vec3{Complex{T}}, 3}: a full vector field in Fourier space;
• pointdata_d::PointData: data associated to vector charges applied on non-uniform points. These are available in pointdata_d.charges and pointdata_d.points;
• ewald_prefactor::Real: the quantity $Γ / V$ where $V$ is the volume of a periodic cell.

The _d suffixes means that data is on the computing device associated to the long-range backend (i.e. on the GPU for GPU-based backends).

Extended help

Implementation details

Fields

All caches must include a common <: LongRangeCacheCommon field which contains common definitions for all backends.

Functions

The following functions must be implemented by a cache:

• transform_to_fourier!,

• interpolate_to_physical!.

NullLongRangeCache <: LongRangeCache

Dummy cache type returned by init_cache_long when long-range computations are disabled.

This is the case when the Ewald splitting parameter $α$ is set to Zero().

LongRangeCacheState

Describes the current state of the long-range fields in Fourier space.

It has two fields, which allow to know which field is currently stored and whether it has been smoothed (Gaussian-filtered) or not:

• quantity::Symbol (can be :undef, :vorticity, :velocity, :streamfunction)

• smoothed::Bool (true if the Gaussian filter associated to Ewald's method has already been applied)

See BiotSavart.get_longrange_field_fourier for more details.

init_cache_long(p::ParamsLongRange, pointdata::PointData, [to::TimerOutput]) -> LongRangeCache

Initialise the cache for the long-range backend defined in p.

Note that, if pc.α === Zero(), then long-range computations are disabled and this returns a NullLongRangeCache.

expected_period(::LongRangeBackend) -> Union{Nothing, Real}

Domain period expected by the backend.

This is used for rescaling input coordinates to the requirements of the backend. For instance, FINUFFT assumes a period $2π$, and therefore coordinates are rescaled if the input data has a period different from $2π$.

folding_limits(::LongRangeBackend) -> Union{Nothing, NTuple{2, Real}}

Domain limits required by the backend.

This is used for folding input coordinates so that they are within the limits expected by the backend. For instance, FINUFFT requires coordinates to be in the $[-3π, 3π]$ interval.

Note that, if a backend defines folding_limits, then it must also define expected_period.

set_num_points!(data::PointData, N::Integer)

Set the total number of non-uniform points that the cache must hold.

This will reallocate space to make all points fit in the cache. It will also reset the contributions of previously-added charges.

add_point_charges!(data::PointData, fs::AbstractVector{<:AbstractFilament}, quad::StaticSizeQuadrature)
add_point_charges!(cache::LongRangeCache, fs::AbstractVector{<:AbstractFilament})

Add vector charges at multiple non-uniform locations.

This can be used for both short-range and long-range computations.

In the case of long-range computations, this must be done before type-1 NUFFTs, to transform from non-uniform data in physical space to uniform data in Fourier space. It must be called before compute_vorticity_fourier!.

has_real_to_complex(::LongRangeBackend) -> Bool
has_real_to_complex(::ParamsLongRange) -> Bool
has_real_to_complex(::LongRangeCacheCommon) -> Bool
has_real_to_complex(::LongRangeCache) -> Bool

Check whether the backend performs real-to-complex (fast) Fourier transforms.

If true, it means that input non-uniform data in physical space can be real-valued, and that uniform data in Fourier space only contains half the total number of modes along the first dimension to account for Hermitian symmetry.

This function is useful in particular for:

• knowing which kind of non-uniform data (vorticities) one must give to the backend;
• knowing how to interpret Fourier-space data, e.g. to compute Fourier spectra.

This function returns false for backends such as FINUFFTBackend, as these require complex input data and don't take advantage of Hermitian symmetry.

compute_vorticity_fourier!(cache::LongRangeCache)

Estimate vorticity in Fourier space.

The vorticity, written to cache.common.uhat_d, is estimated using some variant of non-uniform Fourier transforms (depending on the chosen backend).

Note that this function doesn't perform smoothing over the vorticity using the Ewald operator. Moreover, the resulting vorticity doesn't have the right dimensions, as it must be multiplied by $Γ/V$ (where $Γ$ is the circulation and $V$ is the volume of a periodic cell) to have dimensions $T^{-1}$. In fact, this factor is included in the Ewald operator (see to_smoothed_streamfunction! for details).

Must be called after add_point_charges!.

After calling this function, one may want to use to_smoothed_streamfunction! and/or to_smoothed_velocity! to obtain the respective fields.

to_smoothed_streamfunction!(cache::LongRangeCache)

Convert Fourier-transformed vorticity field to coarse-grained streamfunction field in Fourier space.

This operation can be simply written as:

\[\hat{\bm{ψ}}_{α}(\bm{k}) = ϕ_α(|\bm{k}|) \, \hat{\bm{ω}}(\bm{k})\]

where

\[ϕ_α(k) = \frac{Γ}{V} \, \frac{e^{-k^2 / 4α^2}}{k^2}\]

is the Ewald operator. The effect of this operator is to:

1. invert the Laplacian in $-∇² \bm{ψ} = \bm{\omega}$;
2. smoothen the fields according to the Ewald parameter $α$ (an inverse length scale);
3. rescale values by the vortex circulation $Γ$ and the volume $V$ of a periodic cell so that the streamfunction has the right units ($L^2 T^{-1}$).

This function should be called after compute_vorticity_fourier!. If one only needs the velocity and not the streamfunction, one can also directly call to_smoothed_velocity!.

to_smoothed_velocity!(cache::LongRangeCache)

Convert Fourier-transformed vorticity field to coarse-grained velocity field in Fourier space.

If called right after compute_vorticity_fourier!, this function performs the operation:

\[\hat{\bm{v}}_{α}(\bm{k}) = i \bm{k} × ϕ_α(|\bm{k}|) \, \, \hat{\bm{ω}}(\bm{k}),\]

where $ϕ_α$ is the Ewald operator defined in to_smoothed_streamfunction!.

Optionally, if one is also interested in the streamfunction, one can call to_smoothed_streamfunction! before this function. In that case, the cache already contains the smoothed streamfunction, and only the curl operator ($i \bm{k} ×$) is applied by this function.

interpolate_to_physical!([output::StructVector{<:Vec3},] cache::LongRangeCache)

Perform type-2 NUFFT to interpolate values in cache.common.uhat_d to non-uniform points in physical space.

Results are written to the output vector, which defaults to cache.pointdata_d.charges.

transform_to_fourier!(cache::LongRangeCache)

Transform stored non-uniform data to Fourier space.

This usually corresponds to a type-1 NUFFT.

Non-uniform data must be first added via add_point_charges!.

Base.similar(cache::LongRangeCache, Ns::Dims{3}) -> LongRangeCache

Create a LongRangeCache similar to an existent one, but possibly holding a different amount of wavenumbers in Fourier space.

In principle, the grid resolution Ns can be different from the original one (cache.common.params.Ns). This can be useful for computing high-resolution fields in Fourier space (e.g. for extending the data to higher wavenumbers than allowed by the original cache).

For convenience, point data already in cache.common.pointdata_d is copied to the new cache. This means that, if one already filled the original cache using add_point_charges!, then one can directly call transform_to_fourier! with the new cache to get the vorticity field in Fourier space at the wanted resolution Ns.

ka_generate_kernel(kernel, backend::KA.Backend, x::AbstractArray; [workgroupsize])
ka_generate_kernel(kernel, backend::KA.Backend, ndrange::Dims; [workgroupsize])

Generate statically sized KA kernel.

In this context, "statically sized" means that the kernel will be specifically compiled for the dimensions of the array x, and will be recompiled if an array of a different size is used later.

Here kernel is a KA kernel (a Julia function) annotated with the @kernel macro.

By default, the workgroupsize is determined automatically and may depend on the actual backend (CPU, GPU) and on the array dimensions ndrange = size(x).

KernelAbstractions.get_backend(backend::LongRangeBackend) -> KernelAbstractions.Backend
KernelAbstractions.get_backend(cache::LongRangeCache) -> KernelAbstractions.Backend

Get KernelAbstractions (KA) backend associated to a given long-range backend.

By default this returns KA.CPU(static = true), meaning that things are run on the CPU using threads, and that a static thread assignment is used.

BiotSavartCache

Includes arrays and data required for computation of Biot–Savart integrals.

Fields

• pointdata: contains vector values at points in space. This is used for both short-range and long-range computations;

• shortrange: cache associated to short-range computations;

• longrange: cache associated to long-range computations. It can be NullLongRangeCache() in case the Ewald parameter α was set to Zero();

• to: a TimerOutput instance for measuring the time spent on different functions.

background_vorticity_correction!(
    fields::NamedTuple, fs::AbstractVector{<:AbstractFilament}, params::ParamsBiotSavart,
)

Correct computed fields in case the mean filament vorticity is non-zero.

This ensures that results do not depend on the Ewald splitting parameter $α$ when the filament "charge" is non-zero in a periodic domain.

In practice, as detailed below, this function only corrects the short-range streamfunction, as the velocity is not affected by the background vorticity, and long-range corrections are done implicitly in Fourier space.

Extended help

The mean vorticity associated to the filaments is given by

\[⟨ \bm{ω} ⟩ = -\frac{Γ}{V} ∮_{\mathcal{C}} \mathrm{d}\bm{s},\]

where $V = L^3$ is the periodic domain volume.

This quantity is always zero when all filaments are closed. It can be non-zero if there are infinite unclosed filaments which are not compensated by oppositely-oriented infinite filaments.

When the mean filament vorticity is non-zero, one must compensate it by adding a uniform background vorticity $\bm{ω}_{\text{back}} = -⟨ \bm{ω} ⟩$, so that the total circulation around the volume is zero.

In practice, in the long-range component of Ewald summation, this compensation is done implicitly by setting the zero mode $\hat{\bm{\omega}}(\bm{k} = \bm{0}) = \bm{0}$. As for the short-range component, only the streamfunction needs to be corrected, while the velocity is not affected by the background vorticity.

The correction to the short-range streamfunction is given by

\[\bm{ψ}^{<}_{\text{back}} = \frac{\bm{ω}_{\text{back}}}{4π} ∫_{\mathbb{R}^3} \frac{\operatorname{erfc}(α|\bm{r}|)}{|\bm{r}|} \, \mathrm{d}^3\bm{r} = \frac{\bm{ω}_{\text{back}}}{4 α^2}\]

Therefore, this function simply adds this short-range correction to the streamfunction.