BiotSavart

VortexPasta.BiotSavart Module

BiotSavart

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

source

Biot–Savart parameters

Setting the parameters

VortexPasta.BiotSavart.ParamsBiotSavart Type

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 $|\mathit{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 (1993), 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. See also quadrature_near_singularity for the choice of quadrature rule when this is enabled.

• quadrature_near_singularity = quadrature: quadrature rule to be used when integrating near a singularity when lia_segment_fraction is enabled. By default this is equal to the quadrature rule used for non-local interactions (the quadrature parameter).

Other performance parameters

• use_simd::Bool: whether to use explicit SIMD during the computation of short-range interactions. This can accelerate, in particular, the computation of erfc(αr) appearing in Ewald summation. By default, this is enabled (true) for CPU short-range backends, and there is usually no reason to disable it. On GPU backends (e.g. CellListsBackend(CUDABackend(), …)) this is disabled by default (false). It is not sure whether this could lead to any gains on GPUs, and in fact setting it to true seems to fail on CUDA (but might work on AMD).

• avoid_explicit_erf::Bool = true: whether to avoid explicit computation of erf(αr), which appears when one wants to subtract the local contribution of long-range interactions (the "self"-interaction), and which may be a bit expensive. If this is true, this is avoided by including the local segments in short-range interactions, and then subtracting the full Biot-Savart integrals (instead of long-range ones) evaluated over the local segments. Note that the local integrals (short-range and full) are mathematically singular, and numerically might lead to some loss of accuracy due to catastrophic cancellation, especially when working in single precision (T = Float32).

source

VortexPasta.BiotSavart.autotune Function

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 $\beta$.

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

• inverse splitting parameter $\alpha$;

• cut-off distance $r_{\text{cut}}=\beta /\alpha$;

• Fourier-space cut-off $k_{\text{max}}=2\alpha \beta \approx \pi M/L$, where $M$ is the grid size and $L$ the domain period.

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

Randomness

The parameters selected by this function can be quite random and change from one run to another for a same set of input parameters.

Default accuracy

The default values of β and backend_long correspond to a nominal 6-digit accuracy. See the Extended help for more details.

Periodic domains only

This function only supports periodic domains (finite domain period $L$), since the chosen parameters are irrelevant in the non-periodic case.

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 $\alpha$. 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 $\alpha$, which we write under the form:

$$\alpha =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 $\beta$ and NUFFT parameters

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

Precision digits	$\beta$	NUFFT $w$
3	2.0	2
4	2.5	3
6	3.5	4
8	4.0	5
10	4.5	6
12	5.0	7
14	5.5	8

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 $\sigma =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.

source

Accessing parameters

VortexPasta.BiotSavart.circulation Function

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

Return the circulation Γ associated to each vortex.

source

VortexPasta.BiotSavart.periods Function

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

Return the domain periods in each direction.

source

VortexPasta.BiotSavart.domain_is_periodic Function

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.

source

Exported functions

VortexPasta.BiotSavart.init_cache Function

init_cache(p::ParamsBiotSavart; timer = TimerOutput("BiotSavart")) -> BiotSavartCache

Initialise caches for computing Biot–Savart integrals.

source

VortexPasta.BiotSavart.velocity_on_nodes Function

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

Compute velocity induced by vortex filaments on filament nodes.

Returns a vs vector containing the velocities on filament nodes.

See velocity_on_nodes! for more details.

source

VortexPasta.BiotSavart.velocity_on_nodes! Function

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.

source

VortexPasta.BiotSavart.compute_on_nodes! Function

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(params)
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 is never be called if long-range computations are disabled.

When the callback is called, the available field is simply a Fourier-truncated vorticity, as the coarse-graining due to Ewald's method has not been applied yet.

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)
    (; field, wavenumbers, state,) = BiotSavart.get_longrange_field_fourier(cache)
    @assert state.quantity == :vorticity
    @assert state.smoothing_scale == 0  # unsmoothed field
    # 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
    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)

source

Other convenience functions

VortexPasta.BiotSavart.kelvin_wave_period Function

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

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

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(\lambda )=2\pi /\omega (k)$ where $k=2\pi /\lambda$ is the wavenumber associated to $\lambda$ and $\omega (k)$ is the Kelvin wave dispersion relation:

$$\omega (k)=\frac{\Gamma k^2}{4\pi }[\ln \left(\frac{2}{ka}\right)-\gamma +\frac{1}{2}-\Delta ]$$

where $\gamma \approx 0.5772$ is the Euler–Mascheroni constant.

source

Short-range interactions

Backends

VortexPasta.BiotSavart.ShortRangeBackend Type

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

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

• init_cache_short,

• max_cutoff_distance (optional),

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

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

source

VortexPasta.BiotSavart.NaiveShortRangeBackend Type

NaiveShortRangeBackend <: ShortRangeBackend

Naive computation of short-range interactions.

Maximum cut-off distance

In periodic domains, this backend requires a cut-off distance $r_{\text{cut}}$ not larger than half the domain period $L$ in each direction:

$$r_{\text{cut}}\le \frac{L}{2}$$source

VortexPasta.BiotSavart.CellListsBackend Type

CellListsBackend <: ShortRangeBackend
CellListsBackend([ka_backend = CPU()], [nsubdiv = 1]; [device])

Compute short-range interactions using the cell lists algorithm.

This backend can be significantly faster than the NaiveShortRangeBackend when the cut-off 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 cut-off 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.

Using a GPU

To compute pair interactions on a GPU, pass the corresponding KernelAbstractions.jl backend as the first positional argument.

For example, to use a CUDA device:

using CUDA
backend_short = CellListsBackend(CUDABackend(); kwargs...)

On AMD GPUs the following should work:

using AMDGPU
backend_short = CellListsBackend(ROCBackend(); kwargs...)

If running on a machine with multiple GPU devices, one may use the device keyword argument to choose the device where short-range computations will be performed. This should be a value in 1:ndevices. When using KernelAbstractions.jl, the number of available devices can be obtained using KA.ndevices(ka_backend). If there are 2 or more available devices, then the default is device = 2 (note that device = 1 is used by default for long-range computations, see NonuniformFFTsBackend).

Maximum cut-off distance

The cut-off distance must safisfy the condition:

$$r_{\text{cut}}\le \frac{M}{2M+1}L$$

where $M$ is equal to the nsubdiv parameter. If this is a limitation, one can use the NaiveShortRangeBackend which has a slightly larger limit, $r_{\text{cut}}\le L/2$.

source

Internals

VortexPasta.BiotSavart.ShortRangeCache Type

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;

• pointdata::PointData: updated quadrature locations, charge values and output points (filament nodes);

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

source

VortexPasta.BiotSavart.max_cutoff_distance Function

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.

source

VortexPasta.BiotSavart.init_cache_short Function

init_cache_short(pc::ParamsCommon, p::ParamsShortRange{<:ShortRangeBackend}, pointdata::PointData) -> ShortRangeCache

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

This should be defined for each ShortRangeBackend. In general, implementations should create a copy of pointdata to avoid possible aliasing issues, since the pointdata in the signature is usually the one owned by the full BiotSavartCache.

source

VortexPasta.BiotSavart.add_pair_interactions! Function

add_pair_interactions!(outputs::NamedTuple, cache::ShortRangeCache)

Compute short-range Biot-Savart interactions between pairs of points.

This function evaluates the influence of nearby quadrature points on filament discretisation points. Results are added to the original values in outputs, which means that one should set to zero all values before calling this function (unless one wants to keep previous data, e.g. the velocity associated from other terms not computed here).

The outputs argument is expected to have :velocity and/or :streamfunction fields, containing linear vectors of the same length as the number of output points nodes (cache.pointdata.nodes).

This function does not compute local interactions, i.e. the term associated to local curvature effects in the case of velocity.

source

VortexPasta.BiotSavart.foreach_charge Function

foreach_charge(f::Function, c::ShortRangeCache, x⃗::Vec3)

Apply function f to all charges which are close to point x⃗.

The function should be of the form f(j) where j is the index of a point which is "close" to x⃗.

See also CellLists.foreach_source, which is used when CellListsBackend has been chosen.

source

VortexPasta.BiotSavart.nearby_charges Function

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

Return an iterator over the indices of points that are "close" to the location x⃗.

source

Long-range interactions

Backends

VortexPasta.BiotSavart.LongRangeBackend Type

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;

• 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),

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

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

source

VortexPasta.BiotSavart.NonuniformFFTsBackend Type

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 first positional argument (see below).

Optional arguments

The signature of NonuniformFFTsBackend is:

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

where all arguments besides device 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...)

If running on a machine with multiple GPU devices, one may use the device keyword argument to choose the device where long-range computations will be performed. This should be a value in 1:ndevices. When using KernelAbstractions.jl, the number of available devices can be obtained using KA.ndevices(ka_backend). By default the first device (device = 1) is used for long-range computations.

Keyword arguments

Some relevant keyword arguments which are passed to NonuniformFFTs.jl 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;

• use_atomics = Threads.nthreads() > 4: whether to use atomics (instead of locks) in spreading on the CPU (ignored on GPU devices). See NonuniformFFTs docs for details. By default we enable this for relatively large number of threads, in which case it can be beneficial;

• 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 $\sim {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 (2025):

Precision digits	NUFFT $m$	NUFFT $\sigma$	Ewald $\beta$
3	2	1.5	2.0
4	3	1.5	2.5
6	4	1.5	3.5
8	5	1.5	4.0
10	6	1.5	4.5
12	7	1.5	5.0
14	8	1.5	5.5

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

source

VortexPasta.BiotSavart.ExactSumBackend Type

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.

source

Accessing long-range fields

It may be useful to access computed fields (vorticity, velocity, ...) in Fourier space. For this, one can use the unexported get_longrange_field_fourier function.

VortexPasta.BiotSavart.get_longrange_field_fourier Function

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.

For convenience, if long-range fields were computed on a GPU, this function will also activate the device where these fields are stored (using activate_device!). This is useful when multiple GPUs are being used (e.g. GPU 1 for long-range computations and GPU 2 for short-range ones).

source

Internals

VortexPasta.BiotSavart.PointData Type

PointData{T <: AbstractFloat}
PointData(::Type{T}) -> PointData{T}

Stores point data (values on filaments) used to compute short-range and long-range interactions.

Among the stored fields are:

• nodes::StructVector{Vec3{T}}: filament nodes (discretisation points) where fields like velocity and streamfunction is to be computed;

• points::StructVector{Vec3{T}}: quadrature points on filament segments where "vorticity" is to be evaluated;

• charges::StructVector{Vec3{T}}: "vorticity" vector evaluated on quadrature points. More precisely, this is w * s⃗′ where w is a quadrature weight and s⃗′ is the local tangent vector.

Quadrature points and charges are obtained via interpolation in-between filament nodes.

source

VortexPasta.BiotSavart.LongRangeCache Type

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::NTuple{3, AbstractVector}: Fourier wavenumbers in each direction;

• uhat::StructArray{Vec3{Complex{T}}, 3}: a full vector field in Fourier space;

• pointdata::PointData: data associated to vector charges applied on non-uniform points. These are available in pointdata.charges and pointdata.points;

Note that, for GPU-based backends, practically all arrays (uhat in particular) are GPU arrays, which don't support all operations of CPU arrays (such as for loops).

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!.

source

VortexPasta.BiotSavart.NullLongRangeCache Type

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 $\alpha$ is set to Zero().

source

VortexPasta.BiotSavart.LongRangeCacheState Type

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 its smoothing scale (width $\sigma$ of Gaussian filter):

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

• smoothing_scale::Float64: width $\sigma$ of applied Gaussian filter. This is typically 0 (if the field has not been smoothed), but that can change if one calls to_coarse_grained_vorticity!.

See BiotSavart.get_longrange_field_fourier for more details.

source

VortexPasta.BiotSavart.init_cache_long Function

init_cache_long(p::ParamsBiotSavart, pointdata::PointData) -> LongRangeCache

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

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

source

VortexPasta.BiotSavart.expected_period Function

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, NonuniformFFTs.jl assumes a period $2\pi$, and therefore coordinates are rescaled if the input data has a period different from $2\pi$.

source

VortexPasta.BiotSavart.set_num_points! Function

set_num_points!(data::PointData, Np::Integer, quad::StaticSizeQuadrature)

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.

source

VortexPasta.BiotSavart.add_point_charges! Function

add_point_charges!(data::PointData, fs::AbstractVector{<:AbstractFilament}, params::ParamsBiotSavart)
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 periodic domains, each point is folded into the main periodic cell (in $[0,L{]}^3$). This can be accounted for to speed-up operations done later with point data, in particular in short-range computations when determining the minimal distance between two points in the periodic lattice.

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!.

source

VortexPasta.BiotSavart.has_real_to_complex Function

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.

Currently, this function always returns true. It used to return false for the FINUFFTBackend, which has been removed.

source

VortexPasta.BiotSavart.compute_vorticity_fourier! Function

compute_vorticity_fourier!(cache::LongRangeCache)

Estimate vorticity in Fourier space.

The vorticity, written to cache.common.uhat, is estimated using some kind 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.

Must be called after add_point_charges!.

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

source

VortexPasta.BiotSavart.compute_streamfunction_fourier! Function

compute_streamfunction_fourier!(cache::LongRangeCache)

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

This simply corresponds to inverting a Laplacian ($-{\nabla }^2\mathit{\psi }=\mathit{\omega }$), which in Fourier space is:

$${\hat{\mathit{\psi }}}_{\alpha }(\mathit{k})=\hat{\mathit{\omega }}(\mathit{k})/k²$$

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

source

VortexPasta.BiotSavart.compute_velocity_fourier! Function

compute_velocity_fourier!(cache::LongRangeCache)

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

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

$${\hat{\mathit{v}}}_{\alpha }(\mathit{k})=\frac{i\mathit{k}\times \hat{\mathit{\omega }}(\mathit{k})}{k^2}.$$

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

source

VortexPasta.BiotSavart.to_coarse_grained_vorticity! Function

to_coarse_grained_vorticity!(c::LongRangeCache, ℓ::Real) -> NamedTuple

Compute coarse-grained vorticity in Fourier space.

The vorticity is coarse-grained at scale $\ell$ using a Gaussian filter.

Ideally, for accuracy reasons, the smoothing scale $\ell$ should be larger than the Ewald splitting scale $\sigma =1/\alpha \surd 2$. Also note that applying this function leads to loss of small-scale information, so that one should be careful when performing additional calculations (energy spectra, ...) afterwards.

This can be useful for visualisations or physical analysis. It is not used to compute Biot–Savart velocities, and the scale $\ell$ need not be similar to the splitting lengthscale in Ewald's method.

This function also sets the LongRangeCacheState to quantity = :vorticity and smoothing_scale = ℓ.

This function returns the same thing as get_longrange_field_fourier, including the resulting coarse-grained vorticity in Fourier space.

Example

From a set of vortex filaments, evaluate their self-induced velocity (using compute_on_nodes!) and then the resulting coarse-grained vorticity at their locations:

# First compute the velocity on the filament nodes.
# As an intermediate result, the velocity in Fourier space is stored in the cache.
vs = map(similar ∘ nodes, fs)  # one velocity vector per filament node
fields = (; velocity = vs,)
cache = BiotSavart.init_cache(params)
compute_on_nodes!(fields, cache, fs)

# Now compute the coarse-grained vorticity (result is stored in the cache)
ℓ = 0.1  # smoothing scale
BiotSavart.to_coarse_grained_vorticity!(cache.longrange, ℓ)

# Finally, evaluate the resulting vorticity on filament nodes
ωs_ℓ = map(similar, vs)
BiotSavart.interpolate_to_physical!(cache.longrange)       # interpolate to vortex positions (result stored in the cache)
BiotSavart.copy_long_range_output!(ωs_ℓ, cache.longrange)  # copy results from cache to output array

source

VortexPasta.BiotSavart.interpolate_to_physical! Function

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

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

Results are written to the output vector, which defaults to cache.outputs[1]. This vector is returned by this function.

Using callbacks

The optional callback function should be a function f(û, idx) which takes:

• û::NTuple{3, <:Complex} a Fourier coefficient;

• idx::NTuple{3, Int} the index of û on the Fourier grid.

The function should return a tuple with the same type as û.

This can be used to modify the Fourier-space fields to be interpolated, but without really modifying the values in uhat (and thus the state of the cache). For example, to apply a Gaussian filter before interpolating:

σ = 0.1  # width of Gaussian filter (in physical space)
ks = BiotSavart.get_wavenumbers(cache.longrange)

callback(û::NTuple{3}, idx::NTuple{3}) = let
    k⃗ = @inbounds getindex.(ks, idx)
    k² = sum(abs2, k⃗)
    û .* exp(-σ^2 * k² / 2)
end

interpolate_to_physical!(callback, cache)

GPU usage

As written above, the GPU kernels using the callback might not compile correctly. If that is the case, a possible solution is to use a callable struct instead. For example:

struct GaussianFilter{WaveNumbers} <: Function
    α::Float64
    ks::WaveNumbers
end

@inline function (f::GaussianFilter)(û::NTuple{3}, idx::NTuple{3})
    (; α, ks) = f
    k⃗ = @inbounds getindex.(ks, idx)
    k² = sum(abs2, k⃗)
    û .* exp(-σ^2 * k² / 2)
end

callback = GaussianFilter(0.1, BiotSavart.get_wavenumbers(cache.longrange))  # this is now a callable object
interpolate_to_physical(callback, cache)

If this still doesn't work, one might need to replace broadcasts (.) with calls to map or ntuple.

source

VortexPasta.BiotSavart.transform_to_fourier! Function

transform_to_fourier!(cache::LongRangeCache, prefactor::Real)

Transform stored non-uniform data to Fourier space.

This usually corresponds to a type-1 NUFFT.

The prefactor is a real value that will be used to multiply each non-uniform value (or equivalently each uniform value in Fourier space).

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

source

Base.similar Method

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 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.

source

KernelAbstractions utils

KernelAbstractions.jl is used to write generic code which works on CPUs and different kinds of GPUs.

VortexPasta.BiotSavart.ka_generate_kernel Function

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).

source

KernelAbstractions.get_backend Function

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

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

Note

The word "backend" means two different things here! For KA, it refers to the type of device where kernels are executed (e.g. CPU, CUDABackend, ...).

By default this returns KA.CPU(), meaning that things are run on the CPU using threads.

source

KernelAbstractions.device Function

KernelAbstractions.device(backend::AbstractBackend) -> Int
KernelAbstractions.device(cache::ShortRangeCache) -> Int
KernelAbstractions.device(cache::LongRangeCache) -> Int

Return the device id (in 1:ndevices) where short-range or long-range computations are run.

This can make sense when running on GPUs, where one may want to take advantage of multiple available GPUs on the same machine. On CPUs the device id is generally 1.

source

VortexPasta.BiotSavart.activate_device! Function

BiotSavart.activate_device!(backend::AbstractBackend)
BiotSavart.activate_device!(cache::ShortRangeCache)
BiotSavart.activate_device!(cache::LongRangeCache)

Activate KernelAbstractions device associated to backend or cache.

source

Internals

VortexPasta.BiotSavart.AbstractBackend Type

AbstractBackend

Denotes a "backend" for short-range or long-range Ewald computations.

See ShortRangeBackend and LongRangeBackend for more details.

source

VortexPasta.BiotSavart.BiotSavartCache Type

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.

source

VortexPasta.BiotSavart.background_vorticity_correction! Function

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 $\alpha$ 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.

Note

This function is automatically called by compute_on_nodes! and should never be called manually. This documentation is only included for information purposes.

Extended help

The mean vorticity associated to the filaments is given by

$$⟨\mathit{\omega }⟩=\frac{\Gamma }{V}{\oint }_{\mathcal{C}}\mathrm{d}\mathit{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 ${\mathit{\omega }}_0=-⟨\mathit{\omega }⟩$, so that the total circulation around the volume is zero. This should be interpreted by considering that computations are performed in a rotating reference frame with rotation rate $\mathit{\Omega }=2⟨\mathit{\omega }⟩=-2{\mathit{\omega }}_0$.

In practice, in the long-range component of Ewald summation, this compensation is done implicitly by setting the zero mode $\hat{\mathit{\omega }}(\mathit{k}=\mathbf{0})=\mathbf{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

$${\mathit{\psi }}_0^{<}=\frac{{\mathit{\omega }}_0}{4\pi }{\int }_{{\mathbb{R}}^3}\frac{\mathrm{erfc}(\alpha |\mathit{r}|)}{|\mathit{r}|}{\mathrm{d}}^3\mathit{r}=\frac{{\mathit{\omega }}_0}{4{\alpha }^2}$$

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

source

VortexPasta.BiotSavart.process_point_charges! Function

process_point_charges!(cache::ShortRangeCache)

Process list of point charges.

This is useful for short-range backends like CellListsBackend, which needs to assign a cell to each point charge before finding nearby pairs.

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

source

process_point_charges!(cache::LongRangeCache)

Process list of point charges in long-range cache.

For long-range computations, this should be called after quadrature points and interpolation nodes have been set, either via add_point_charges!(cache, fs) or by directly modifying cache.pointdata. It must be called before any calls to compute_vorticity_fourier! or interpolate_to_physical!.

This function will process (and possibly modify) data in cache.pointdata. Therefore, it should only be called once after points have been set. For example, the NonuniformFFTsBackend assumes a domain of size [0, 2π]ᵈ, and thus a point transformation is needed if the domain has a different size.

source

VortexPasta.BiotSavart.copy_output_values_on_nodes! Function

copy_output_values_on_nodes!([op::Function], vs::AbstractVector, vs_d::AbstractVector, [vs_h::AbstractVector])

Copy computed values onto a vector of vectors.

Parameters

• op: an optional binary operation op(old, new) to be applied when writing each value into vs. An usual choice is +, which means that copied values will be added to previously existing values. By default this is op(old, new) = new, meaning that old values are replaced with new values.

• vs: output vector, usually a vector of vectors used to store velocity or streamfunction on filament nodes.

• vs_d: vector where quantities (e.g. velocity) where computed. Unlike vs, this is a "linear" vector, meaning that there is a single vector of quantities (Vec3) for all filament nodes across all filaments. Note that this vector may also be on a GPU device (_d = device).

• vs_h: this is an optional CPU vector (_h = host), which may be passed if vs_d is not a CPU array. In that case, it will be used as an intermediate buffer array, and thus passing it may help reduce CPU allocations.

source