Diagnostics

VortexPasta.Diagnostics Module

Diagnostics

Contains tools for computing different diagnostics (total energy, energy spectra, ...) from simulation data.

source

Kinetic energy

VortexPasta.Diagnostics.kinetic_energy Function

kinetic_energy(iter::VortexFilamentSolver; quad = nothing) -> Real
kinetic_energy(fs, ψs, Γ::Real, [Ls]; quad = nothing) -> Real
kinetic_energy(fs, ψs, p::ParamsBiotSavart; quad = nothing) -> Real

Compute kinetic energy of velocity field induced by a set of vortex filaments.

This function simply calls kinetic_energy_from_streamfunction.

source

VortexPasta.Diagnostics.kinetic_energy_from_streamfunction Function

kinetic_energy_from_streamfunction(iter::VortexFilamentSolver; quad = nothing)
kinetic_energy_from_streamfunction(fs, ψs, Γ::Real, [Ls]; quad = nothing)
kinetic_energy_from_streamfunction(fs, ψs, p::ParamsBiotSavart; quad = nothing)

Compute kinetic energy per unit mass (units $L^2{T}^{-2}$) from streamfunction values at filament nodes in a periodic domain.

The kinetic energy per unit mass of the velocity field induced by a set of vortex filaments can be obtained as:

$$E=\frac{\Gamma }{2V}\oint \mathit{\psi }(\mathit{s})\cdot \mathrm{d}\mathit{s}$$

where $\Gamma$ is the vortex circulation and $V$ is the volume of interest (e.g. the volume of a periodic cell in periodic domains).

Arguments

Mandatory arguments

• fs: vortex filament locations;

• ψs: streamfunction values at filament nodes;

• Γ::Real: quantum of circulation; or

• p::ParamsBiotSavart: Biot–Savart parameters (see ParamsBiotSavart).

Optional arguments

• Ls::Tuple = (Lx, Ly, Lz): domain size in each direction. If not given, the domain volume $V$ is taken to be 1 (see Non-periodic domains below).

Optional keyword arguments

• quad = nothing: optional quadrature rule (e.g. quad = GaussLegendre(4)) used to evaluate line integrals. If nothing, only values at nodes are used (cheaper). Otherwise, if a quadrature rule is passed, interpolations are performed and extra allocations are needed.

Non-periodic domains

If the domain is not periodic, that is, if one or more values in Ls is Infinity, then the domain volume $V$ in the expression above will be set to 1. This corresponds to computing an energy per unit density (units $L^5{T}^{-2}$) instead of per unit mass, meaning that one should multiply by the fluid density $\rho$ ($M{L}^{-3}$) to get an actual kinetic energy ($M{L}^2{T}^{-2}$).

Note that in non-periodic domains one may also use the kinetic_energy_nonperiodic function, which uses a different definition commonly used for open domains and which does not require streamfunction values (but only velocity values on filament nodes). However, energy computed using that definition may not be properly conserved when it should.

Therefore, it is recommended to always use kinetic_energy_from_streamfunction, even in non-periodic domains.

source

VortexPasta.Diagnostics.kinetic_energy_nonperiodic Function

kinetic_energy_nonperiodic(iter::VortexFilamentSolver; quad = nothing) -> Real
kinetic_energy_nonperiodic(fs, vs, Γ::Real; quad = nothing) -> Real
kinetic_energy_nonperiodic(fs, vs, p::ParamsBiotSavart; quad = nothing) -> Real

Compute kinetic energy per unit density (units $L^5{T}^{-2}$) from velocity values at filament nodes in an open (non-periodic) domain.

This function returns the kinetic energy over the infinite fluid volume, which only makes sense in an open domain such that the velocity tends to zero far from the vortices.

In an open domain, the kinetic energy of the velocity field induced by a set of vortex filaments can be obtained as:

$$E=\rho \Gamma \oint \mathit{v}\cdot (\mathit{s}\times \mathrm{d}\mathit{s})$$

where $\Gamma$ is the vortex circulation and $\rho$ is the fluid density. This definition assumes that the velocity field tends to zero far from the vortices. This is true in open domains, but not in periodic ones.

Energy conservation

Energy computed using this definition may present small temporal fluctuations in cases where energy should be conserved. For this reason, it is recommended to always use kinetic_energy_from_streamfunction, which displays proper energy conservation properties. In general, this definition will slightly overestimate the energy obtained from the streamfunction.

This function returns the energy per unit density $E/\rho$.

Mandatory arguments

• vs: velocity values at filament nodes;

• fs: vortex filament locations;

• Γ::Real: quantum of circulation.

Optional keyword arguments

See kinetic_energy_from_streamfunction for details.

Periodic domains

In periodic domains this function will give wrong results, since there are boundary terms coming from integration by parts which are neglected in the above definition (assuming the velocity goes to zero far from the vortices, which is not the case in a periodic domain).

In this case, the kinetic_energy_from_streamfunction function should be used instead.

source

Energy injection rate

VortexPasta.Diagnostics.energy_injection_rate Function

energy_injection_rate(iter::VortexFilamentSolver, [vL]; quad = nothing) -> Real
energy_injection_rate(fs, vL, vs, p::ParamsBiotSavart; quad = nothing) -> Real

Compute energy injection rate from current filament velocities.

This is typically the energy injected by an advecting vortex velocity vL, given the self-induced (Biot–Savart) velocity vs. It can be negative if vL actually dissipates energy. It does not include the contributions to the dissipation by vortex reconnections or numerical resolution effects.

Arguments

• fs: vortex filament locations;

• vL: velocity of filament nodes (including forcing or external velocities);

• vs: self-induced filament velocity (from Biot–Savart law);

• p::ParamsBiotSavart: Biot–Savart parameters (see ParamsBiotSavart);

• quad: an optional quadrature rule (such as GaussLegendre) for accurate line integrals.

Extended help

Definition

The energy injection rate (in $L^2{T}^{-3}$ units) can be expressed as

$${\epsilon }_{\text{inj}}=\frac{\Gamma }{V}\oint \left[{\mathit{s}}^{\prime }\times {\mathit{v}}_{\text{s}}\right]\cdot {\mathit{v}}_{\text{L}}\mathrm{d}\xi ,$$

where ${\mathit{v}}_{\text{s}}$ is the self-induced vortex velocity (according to the Biot–Savart law) and ${\mathit{v}}_{\text{L}}$ is the actual velocity used to advect the vortices. Typically, the latter can be written as

$${\mathit{v}}_{\text{L}}={\mathit{v}}_{\text{s}}+{\mathit{v}}_{\text{ext}}$$

where ${\mathit{v}}_{\text{ext}}$ is an externally applied velocity (for example, representing the interaction with a normal fluid). Note that ${\mathit{v}}_{\text{s}}$ doesn't contribute to energy injection, so that energy is conserved (excluding other dissipative effects) in the absence of an external velocity.

Energy transfer from "singular" to total velocity

It is also possible to estimate an energy transfer from the "singular" part of the BS velocity to its non-singular part. Here we define the "singular" BS velocity as ${\mathit{v}}_{\text{s}}^{\infty }\equiv \frac{\Gamma }{4\pi }{\mathit{s}}^{\prime }\times {\mathit{s}}^{″}$. Choosing this as ${\mathit{v}}_{\text{L}}$ leads to:

$${\epsilon }_{\text{inj}}^{\infty }=\frac{{\Gamma }^2}{4\pi V}\oint {\mathit{s}}^{″}\times {\mathit{v}}_{\text{s}}\mathrm{d}\xi .$$

With a negative sign, this may be interpreted as the energy transfered from the non-local to the local part of the Biot–Savart velocity. Note that this is directly related to the filament stretching rate (see stretching_rate).

To compute this estimate, one should pass vL = CurvatureVector() as the second argument to this function.

source

Helicity

VortexPasta.Diagnostics.helicity Function

helicity(iter::VortexFilamentSolver; quad = nothing) -> Real
helicity(fs, vs, Γ::Real; quad = nothing) -> Real
helicity(fs, vs, p::ParamsBiotSavart; quad = nothing) -> Real

Compute helicity of a vortex configuration.

The helicity is defined as:

$$\mathcal{H}=\int \mathit{u}(\mathit{x})\cdot \mathit{\omega }(\mathit{x}){\mathrm{d}}^3\mathit{x}=\Gamma \oint \mathit{v}(\mathit{s})\cdot \mathrm{d}\mathit{s}$$

where $\mathit{u}(\mathit{x})$ and $\mathit{\omega }(\mathit{x})=\nabla \times \mathit{u}(\mathit{x})$ are the velocity and vorticity fields, and $\mathit{v}(\mathit{s})$ is the velocity of the set of vortex filaments.

Note that the returned helicity has units of a squared circulation (${\Gamma }^2\sim L^4{T}^{-2}$).

Arguments

• fs: list of vortex filaments (or a single AbstractFilament);

• vs: list of velocities on filament nodes (or a single vector of velocities on a filament);

• p: Biot–Savart parameters; or

• Γ: vortex circulation.

See e.g. kinetic_energy_from_streamfunction for the meaning of the optional quad keyword argument.

source

Filament length

See filament_length in the Filaments module. For convenience, filament_length is re-exported by Diagnostics, meaning that one can do:

using VortexPasta.Diagnostics
filament_length(...)

without needing to import VortexPasta.Filaments.

Stretching rate

VortexPasta.Diagnostics.stretching_rate Function

stretching_rate(iter::VortexFilamentSolver; quad = nothing) -> Real
stretching_rate(fs, vs; quad = nothing) -> Real

Compute stretching rate of one or more vortices.

The stretching rate has units $L{T}^{-1}$ and is given by:

$$\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}t}=\oint \frac{\partial \mathit{v}}{\partial \xi }\cdot \mathrm{d}\mathit{s}=-\oint \mathit{v}\cdot {\mathit{s}}^{″}\mathrm{d}\xi$$

where $\xi$ is the arc length, $\mathit{v}(\xi )$ the local filament velocity, ${\mathit{s}}^{″}(\xi )$ the local curvature vector, and $\mathcal{L}$ the instantaneous vortex length. The last equality is obtained using integration by parts.

In the implementation, the last expression is the one used to compute the stretching rate.

Mandatory arguments

• vs: velocity values at filament nodes;

• fs: vortex filament locations.

Optional keyword arguments

• quad = nothing: optional quadrature rule (e.g. quad = GaussLegendre(4)) used to evaluate line integrals. If nothing, only values at nodes are used (cheaper).

source

Vortex impulse

VortexPasta.Diagnostics.vortex_impulse Function

vortex_impulse(iter::VortexFilamentSolver; quad = nothing) -> Vec3
vortex_impulse(f; quad = nothing) -> Vec3

Estimate normalised impulse of one or more vortex filaments.

The vortex impulse is defined as

$$\mathit{p}=\frac{1}{2}\rho \Gamma \oint \mathit{s}\times \mathrm{d}\mathit{s}$$

where $\rho$ is the fluid density and $\Gamma$ the circulation about the vortex. Note that this function returns the normalised impulse $\mathit{p}/\rho \Gamma$. The returned impulse has units of $L^2$ (an area).

Note that, for a circular vortex ring of radius $R$, its impulse is $\mathit{p}=\rho \Gamma A$ where $A=\pi R^2$ is the area enclosed by the ring (and the orientation is equal to its direction of propagation, i.e. normal to the plane where the ring lives).

source

Energy and helicity spectrum

VortexPasta.Diagnostics.init_spectrum Function

Diagnostics.init_spectrum(iter::VortexFilamentSolver) -> (ks, Ek)
Diagnostics.init_spectrum(cache) -> (ks, Ek)

Initialise fields for storing an energy or helicity spectrum.

Returns a wavenumber vector ks and an uninitialised spectrum Ek with the right dimensions, which can be then passed to energy_spectrum! or `helicity_spectrum!.

The returned arrays are always on the CPU, even when the cache contains GPU data.

See energy_spectrum! for details on the cache argument.

source

VortexPasta.Diagnostics.energy_spectrum Function

energy_spectrum(iter::VortexFilamentSolver) -> (ks, Ek)
energy_spectrum(cache) -> (ks, Ek)

Compute kinetic energy spectrum associated to vortex filament state.

Returns a tuple of vectors (ks, Ek) where ks contains the probed wavenumbers and Ek the energy associated to each wavenumber.

See also energy_spectrum! for a non-allocating variant and for more details.

source

VortexPasta.Diagnostics.energy_spectrum! Function

energy_spectrum!(Ek::AbstractVector, ks::AbstractVector, cache)

Compute kinetic energy spectrum associated to vortex filament state.

Here cache contains the results of long-range Biot–Savart computations. It can be either:

• a LongRangeCache;

• a BiotSavartCache (which contains a LongRangeCache);

• a VortexFilamentSolver from the Timestepping module (which contains a BiotSavartCache).

The energy spectrum is computed from a recent Biot–Savart calculation using fast Ewald summation. More precisely, it is computed from the (Fourier-truncated) velocity field in Fourier space. The LongRangeCache associated to the calculation is expected to currently contain this field.

The cache can also contain the vorticity field in Fourier space (the result obtained right after performing a NUFFT from the filament locations, see BiotSavart.compute_vorticity_fourier!. In this case this function does the right thing and also computes the spectrum of the associated velocity field.

The vectors Ek and ks are expected to have the same length. Moreover, the vector of wavenumbers ks should satisfy ks[begin] == 0 and have a constant step Δk = ks[i + 1] - ks[i]. For convenience, the init_spectrum function can be used to create these vectors.

See also energy_spectrum for an allocating variant which doesn't need predefined Ek and ks vectors.

source

VortexPasta.Diagnostics.helicity_spectrum Function

helicity_spectrum(iter::VortexFilamentSolver) -> (ks, Hk)
helicity_spectrum(cache) -> (ks, Hk)

Compute helicity spectrum associated to vortex filament state.

Returns a tuple of vectors (ks, Hk) where ks contains the probed wavenumbers and Hk the helicity associated to each wavenumber.

See also helicity_spectrum! for a non-allocating variant and for more details.

source

VortexPasta.Diagnostics.helicity_spectrum! Function

helicity_spectrum!(Hk::AbstractVector, ks::AbstractVector, cache)

Compute helicity spectrum associated to vortex filament state.

For a single wave vector $\mathit{k}$, the associated helicity is defined as

$$H(\mathit{k})\equiv \frac{V}{\Delta k}{\mathit{v}}_{\mathit{k}}^{\ast }\cdot {\mathit{\omega }}_{\mathit{k}}=2\frac{V}{\Delta k}\mathit{k}\cdot \left[{\mathit{v}}_{\text{r}}\times {\mathit{v}}_{\text{i}}\right].$$

In the last expression, the subscripts $\text{r}$ and $\text{i}$ respectively denote real and imaginary parts. That result can be easily found by doing some algebra, after replacing ${\mathit{\omega }}_{\mathit{k}}=\mathrm{i}\mathit{k}\times {\mathit{u}}_{\mathit{k}}$. It explicitly shows that $H(\mathit{k})$ is purely real.

As usual, this function computes the sum of $H(\mathit{k})$ over spherical shells of radius $k_n\le |\mathit{k}|<k_n+\Delta k=k_{n+1}$.

Defined as above, the total helicity in a (cubic) periodic domain of volume $V=L^3$ is $H={\sum }_{\mathit{k}}H(\mathit{k})\Delta k$ where $\Delta k=2\pi /L$ is the resolution in wavenumber space.

See energy_spectrum! for details on accepted cache arguments.

source

Spectral energy fluxes

VortexPasta.Diagnostics.energy_flux Function

energy_flux(iter::VortexFilamentSolver, Nk::Integer; quad = nothing) -> (ks, fluxes)
energy_flux(iter::VortexFilamentSolver, ks::AbstractVector; quad = nothing) -> (ks, fluxes)

Compute the energy "flux" across scales due to various velocity terms.

The second argument can be either the wanted number of output wavenumbers Nk, or the actual output wavenumbers ks. In the first case, this function will choose logarithmically-spaced wavenumbers ks across which the fluxes will be computed. In general, the length of ks will be slightly lower than Nk to avoid possible repetitions.

Note that the cost of computing the fluxes across each wavenumber is very high (we basically perform one NUFFT per wavenumber), so Nk should not be too large, and this function should not be called too often.

This function returns:

• ks::AbstractVector: a vector of wavenumbers;

• fluxes::NamedTuple: a list of energy fluxes across scales due to different velocity terms. Each element is a vector of real values. For example, fluxes.vs is a vector containing the energy flux across scales due to the Biot–Savart velocity (see Extended help below).

Extended help

Depending on the fields available in iter, the returned fluxes can include the fields:

• vs: energy flux across scale $k$ due to the Biot–Savart self-induced velocity;

• vinf: estimated non-local energy flux from scales $\le k$ to scales $k\sim \infty$;

• vf: accumulated energy injection rate (by external forcing term) up to scale $k$;

• vdiss: accumulated energy dissipation rate (by external dissipation term) up to scale $k$.

Note that all of these fluxes are properly signed according to the usual convention (i.e. they are ready to plot). In particular, both vf and vdiss are positive if the respective terms inject and dissipate energy as expected, while vs is positive in the presence of a direct energy cascade.

Definitions

The energy flux (fluxes.vs) across scale $k$ is defined as:

$${\Pi }_k=-\frac{\Gamma }{V}\oint \left[{\mathit{s}}^{\prime }\times {\mathit{v}}_{\text{s}}^{<}\right]\cdot {\mathit{v}}_{\text{s}}\mathrm{d}\xi ,$$

where ${\mathit{v}}_{\text{s}}$ is the self-induced vortex velocity (according to the Biot–Savart law) and the superscript $<$ denotes a low-pass filter up to scale $k$. Note that this looks very similar to the energy_injection_rate definition, and this is because ${\Pi }_k$ can indeed be interpreted as the dissipation of large-scale energy (associated to ${\mathit{v}}_{\text{s}}^{<}$) by the full Biot–Savart velocity ${\mathit{v}}_{\text{s}}$.

Similarly, the estimated non-local energy flux (fluxes.vinf) from scales $\le k$ to scales $k\sim \infty$ is defined as:

$$\begin{array}{rl}{\Pi }_k^{\infty }& =-\frac{\Gamma }{V}\oint \left[{\mathit{s}}^{\prime }\times {\mathit{v}}_{\text{s}}^{<}\right]\cdot {\mathit{v}}_{\text{s}}^{\infty }\mathrm{d}\xi \\ & =-\frac{{\Gamma }^2}{4\pi V}\oint {\mathit{s}}^{″}\cdot {\mathit{v}}_{\text{s}}^{<}\mathrm{d}\xi ,\end{array}$$

where ${\mathit{v}}_{\text{s}}^{\infty }\equiv \frac{\Gamma }{4\pi }{\mathit{s}}^{\prime }\times {\mathit{s}}^{″}$ is an estimate for the "singular" part of the self-induced velocity, associated to very high wavenumbers $k$. This is only an estimate which allows to have an idea of the non-local energy transfers from large to small scales. Also note the close similarity of this expression with the filament stretching rate (see stretching_rate). Indeed, ${\Pi }_k^{\infty }$ is proportional to the increase of filament length due to the large-scale velocity ${\mathit{v}}_{\text{s}}^{<}$.

If an external dissipation term is included, the accumulated energy dissipation rate (fluxes.vdiss) is given by:

$$D_k=-\frac{\Gamma }{V}\oint \left[{\mathit{s}}^{\prime }\times {\mathit{v}}_{\text{s}}^{<}\right]\cdot {\mathit{v}}_{\text{diss}}\mathrm{d}\xi ,$$

where ${\mathit{v}}_{\text{diss}}$ is the dissipative velocity.

Finally, the accumulated energy injection rate (fluxes.vf) is (note the $+$ sign!):

$$F_k=+\frac{\Gamma }{V}\oint \left[{\mathit{s}}^{\prime }\times {\mathit{v}}_{\text{s}}^{<}\right]\cdot {\mathit{v}}_{\text{f}}\mathrm{d}\xi ,$$

where ${\mathit{v}}_{\text{f}}$ is the forcing velocity.

source

VortexPasta.Diagnostics.energy_transfer_matrix Function

energy_transfer_matrix(iter::VortexFilamentSolver, Nk::Integer; quad = nothing) -> (ks, transfers)
energy_transfer_matrix(iter::VortexFilamentSolver, ks::AbstractVector; quad = nothing) -> (ks, transfers)

Compute energy transfers between different Fourier shells.

Returns a square antisymmetric matrix of dimensions (Nk + 1, Nk + 1) where Nk is the length of the returned ks vector (which might be smaller than the input Nk, if there are not enough resolved wavenumbers). Note that the ks vector does not include the mode k = 0.

This function defines Nk + 1 shells, such that the i-th shell contains wavenumbers ks[i - 1] < k ≤ ks[i], with the conventions ks[0] = 0 and ks[Nk + 1] = ∞. In other words, the last "shell" contains all small wavenumbers k > ks[end].

See energy_flux for more details on the accepted arguments.

Definitions

The energy transfer between Fourier shells A and B is defined as

$$T_{AB}=\frac{\Gamma }{V}\oint \left[{\mathit{s}}^{\prime }\times {\mathit{v}}_{\text{s}}^{\text{A}}\right]\cdot {\mathit{v}}_{\text{s}}^{\text{B}}\mathrm{d}\xi ,$$

where ${\mathit{v}}_{\text{s}}^{\text{A}}$ and ${\mathit{v}}_{\text{s}}^{\text{B}}$ are Fourier band-pass-filtered fields, including only the coefficients within the chosen shells.

source

Length scales

VortexPasta.Diagnostics.integral_lengthscale Function

integral_lengthscale(ks::AbstractVector, Ek::AbstractVector, Etot::Real, Lvort::Real, p::ParamsBiotSavart) -> Real
integral_lengthscale(ks::AbstractVector, Ek::AbstractVector, Etot::Real, Lvort::Real, iter::VortexFilamentSolver) -> Real
integral_lengthscale(ks::AbstractVector, Ek::AbstractVector, [Etot::Real]) -> Real

Integral length scale associated to velocity field.

The first two variants (taking Etot and Lvort parameters) should be preferred to compute the integral length scale of a superfluid flow.

The last one may be used for smooth velocity fields such that the energy spectrum quickly decays to zero. Note that, if the total kinetic energy Etot is not passed, then it is directly estimated from the spectrum.

Dependencies

Before calling this function, one first needs to:

• call Diagnostics.energy_spectrum (or Diagnostics.energy_spectrum!) to obtain ks and Ek;

• call Diagnostics.kinetic_energy to obtain Etot (total kinetic energy);

• call Diagnostics.filament_length to obtain Lvort (total vortex length).

Definition

The integral length scale can be defined from the energy spectrum $E(k)$ as

$$L_0=\frac{\pi }{2v_{\text{rms}}^2}{\int }_0^{\infty }\frac{E(k)}{k}\mathrm{d}k=\frac{3\pi }{4E}{\int }_0^{\infty }\frac{E(k)}{k}\mathrm{d}k$$

(see e.g. Pope 2000, eq. 6.225). Note that the velocity variance can be related to the kinetic energy $E$ (Etot) simply by $E=(3/2)v_{\text{rms}}^2$ in three dimensions.

For a non-smooth velocity field such as the one generated by vortex lines, $E(k)$ decays slowly (as $k^{-1}$). Therefore, if we truncate the above integral (which converges) at some maximum wavenumber $k_{\text{max}}$, we necessarily incur some truncation error which may be non-negligible. For this reason, we estimate the contribution beyond $k_{\text{max}}$ by assuming that the spectrum takes the form $E(k)=({\Gamma }^2{L}_v/4\pi V)k^{-1}$, where $L_v$ is the total vortex length (Lvort) and $V$ the volume of the periodic domain (see e.g. Polanco (2025) supplementary material).

source