Timestepping

Timestepping

Module defining timestepping solvers for vortex filament simulations.

VortexFilamentProblem(fs::AbstractVector{<:AbstractFilament}, tspan::NTuple{2}, p::ParamsBiotSavart)

Define a vortex filament problem.

Arguments:

• fs: initial vortex positions;

• tspan = (t_begin, t_end): time span;

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

See init for initialising a solver from a VortexFilamentProblem.

init(prob::VortexFilamentProblem, scheme::TemporalScheme; dt::Real, kws...) -> VortexFilamentSolver

Initialise vortex filament problem.

Returns a VortexFilamentSolver which can be advanced in time using either step! or solve!.

Mandatory positional arguments

• prob: a VortexFilamentProblem containing the problem definition.

• scheme: a timestepping scheme (see TemporalScheme).

Mandatory keyword arguments

• dt::Real: the simulation timestep.

Optional keyword arguments

• alias_u0 = false: if true, the solver is allowed to modify the initial vortex filaments.

• refinement = NoRefinement(): criterion used for adaptive refinement of vortex filaments. See RefinementCriterion for a list of possible criteria.

• reconnect = NoReconnections(): criterion used to perform vortex reconnections. See ReconnectionCriterion for a list of possible criteria.

• adaptivity = NoAdaptivity(): criterion used for adaptively setting the timestep dt. See AdaptivityCriterion for a list of possible criteria.

• dtmin = 0.0: minimum dt for adaptive timestepping. If dt < dtmin, the solver is stopped with an error.

• fast_term = LocalTerm(): for IMEX and multirate schemes, this determines what is meant by "fast" and "slow" dynamics. This can either be LocalTerm or ShortRangeTerm.

• LIA = false: if true, only use the local induction approximation (LIA) to advance vortex filaments, ignoring all non-local interactions. Note that reconnections may still be enabled.

• fold_periodic = true: if true (default), vortices will be recentred onto the main unit cell when using periodic boundary conditions. It may be convenient to disable this for visualisation purposes. This setting doesn't affect the results (velocity of each filament, reconnections, …), besides possible spatial translations of the filaments proportional to the domain period.

• callback: a function to be called at the end of each timestep. The function must accept a single argument iter::VortexFilamentSolver. Filaments should not be modified by this function. In that case use affect! instead. See notes below for more details.

• affect!: similar to callback, but allows to modify filament definitions. See notes below for more details.

• timer = TimerOutput("VortexFilament"): an optional TimerOutput for recording the time spent on different functions.

• external_velocity / external_streamfunction: allows to add an external velocity field to "force" the filaments. See "Adding an external velocity" below for more details.

• stretching_velocity: allows to add an external "stretching" velocity to the filaments. This can be considered as an external energy (or length) injection mechanism. However, it can also lead to spurious Kelvin waves. See "Adding a stretching velocity" below for more details.

• forcing: allows to modify the vortex velocities due to mutual friction with an imposed normal fluid velocity field. See "Forcing via a normal fluid" below for more details.

Extended help

Difference between callback and affect!

The difference between the callback(iter) and affect!(iter) functions is the time at which they are called:

• the affect! function is called before performing Biot-Savart computations from the latest filament positions. In other words, the fields in iter.quantities are not synchronised with iter.fs, and it generally makes no sense to access them. Things like energy estimates will be incorrect if done in affect!. On the other hand, the affect! function allows to modify iter.fs before Biot–Savart computations are performed. In particular, to inject new filaments, call inject_filament! from within the affect! function.

• the callback function is called after performing Biot-Savart computations. This means that the fields in iter.quantities have been computed from the latest filament positions. However, one must not modify iter.fs, or otherwise the velocity iter.vs (which will be used at the next timestep) will no longer correspond to the filament positions.

Forcing / energy injection methods

There are multiple implemented methods for affecting and injecting energy onto the vortex system. These methods can represent different physical mechanisms (or in fact be quite unphysical), and some of them can make more sense than others when simulating e.g. zero-temperature superfluids.

Adding an external velocity

One can set the external_velocity keyword argument to impose an external velocity field $\bm{v}_{\text{f}}(\bm{x}, t)$. In this case, the total velocity at a point $\bm{x}$ (which can be on a vortex filament) will be

\[\bm{v}(\bm{x}) = \bm{v}_{\text{BS}}(\bm{x}) + \bm{v}_{\text{f}}(\bm{x})\]

where $\bm{v}_{\text{BS}}$ is the velocity obtained from the Biot–Savart law.

This can be a way of applying an external forcing (and injecting energy) to the vortex system.

The external velocity should be given as a function, which should have the signature v_ext(x⃗::Vec3, t::Real) -> Vec3. The function should return a Vec3 with the velocity at the location x⃗ and time t. For example, to impose a constant sinusoidal forcing in the $x$ direction which varies along $z$, the external_velocity keyword argument could look like:

external_velocity = (x⃗, t) -> Vec3(0.1 * sin(2 * x⃗.z), 0, 0)

One usually wants the external velocity to be divergence-free, to preserve the incompressibility of the flow.

Notes on energy estimation

When setting an external velocity $\bm{v}_{\text{f}}(\bm{x}, t)$, one may also want to impose an external streamfunction field $\bm{ψ}_{\text{f}}(\bm{x}, t)$ via the external_streamfunction keyword argument. Such a field plays no role in the dynamics, but it allows to include the kinetic energy associated to the external velocity field when calling Diagnostics.kinetic_energy_from_streamfunction. If provided, the external streamfunction should satisfy $\bm{v}_{\text{f}} = \bm{∇} × \bm{ψ}_{\text{f}}$.

More precisely, when applying an external velocity field, the total kinetic energy is

\[\begin{align*} E &= \frac{1}{2V} ∫ |\bm{v}_{\text{BS}} + \bm{v}_{\text{f}}|^2 \, \mathrm{d}^3\bm{x} \\ &= \frac{1}{2V} \left[ Γ ∮_{\mathcal{C}} \left( \bm{ψ}_{\text{BS}} + 2 \bm{ψ}_{\text{f}} \right) ⋅ \mathrm{d}\bm{s} + ∫ |\bm{v}_{\text{f}}|^2 \, \mathrm{d}^3\bm{x} \right], \end{align*}\]

where $\bm{ψ}_{\text{BS}}$ is the streamfunction obtained from the Biot-Savart law. Accordingly, when passing external_streamfunction, the streamfunction values attached to vortex filaments (in iter.ψs) will actually contain $\bm{ψ}_{\text{BS}} + 2 \bm{ψ}_{\text{f}}$ (note the factor 2) to enable the computation of the full kinetic energy. Moreover, calling Diagnostics.kinetic_energy_from_streamfunction with a VortexFilamentSolver will automatically include an estimation of the kinetic energy of the external velocity field (the volume integral above).

Adding a stretching velocity

One can set the stretching_velocity keyword argument to impose a stretching velocity $\bm{v}_\text{L}(ξ, t)$ on filament locations $\bm{s}(ξ, t)$. This velocity will be parallel (and generally opposite) to the curvature vector $\bm{s}'' = ρ \bm{N}$, where $ρ$ is the curvature and $\bm{N}$ the unit normal.

More precisely, the stretching_velocity parameter allows to set a filament velocity of the form

\[\bm{v}_{\text{L}}(ξ) = -v_{\text{L}}[ρ(ξ)] \, \bm{N}(ξ)\]

where $v_{\text{L}}$ is a velocity magnitude which can optionally depend on the local curvature value $ρ$. The stretching_velocity parameter should therefore be a function v(ρ::Real) -> Real.

Note that $\bm{v} ⋅ \bm{s}''$ can be roughly interpreted as a local stretching rate (this is actually true for the integrated quantity, see Diagnostics.stretching_rate). Therefore, one may want $v_{\text{L}}$ to be inversely proportional to the local curvature $ρ$, so that the local stretching rate is approximately independent of the filament location $ξ$. To achieve this, one may set

stretching_velocity = ρ -> min(γ / ρ, v_max)

where γ is a constant (units $T^{-1}$) setting the stretching magnitude, and v_max is a maximum stretching velocity used to avoid very large velocities in low-curvature regions. One could do something fancier and replace the min with a smooth regularisation such as

stretching_velocity = ρ -> -expm1(-ρ / ρ₀) * (γ / ρ)  # note: expm1(x) = exp(x) - 1

for some small curvature $ρ₀$. The maximum allowed velocity will then be vmax = γ / ρ₀.

Forcing via a normal fluid

Alternatively to the above approaches, one may impose a "normal" fluid velocity field to modify, via a mutual friction force, the velocity of the vortex filaments. To do this, one first needs to construct a NormalFluidForcing (see that link for definitions and examples) representing a normal fluid velocity field. Then, the obtained object should be passed as the optional forcing argument.

In principle, this method can be combined with the previous ones. The normal fluid forcing will be applied after all the other forcing methods. In other words, the vortex velocity $\bm{v}_{\text{s}}$ used to compute the mutual friction velocity includes all the other contributions (i.e. external and stretching velocities).

There is also a FourierBandForcing type which behaves similarly to NormalFluidForcing but is more localised in scale space, as it uses a band-pass filtered version of the superfluid velocity induced by the vortices.

Injecting filaments over time

Another way of injecting energy is simply by adding vortices to the simulation from time to time. This can be achieved by using an affect! function. See inject_filament! for some more details.

VortexFilamentSolver

Contains the instantaneous state of a vortex filament simulation.

Must be constructed using init.

Included fields

Some useful fields included in a VortexFilamentSolver are:

• prob: associated VortexFilamentProblem (including Biot–Savart parameters);

• fs: current state of vortex filaments in the system;

• quantities: a NamedTuple containing data on filament nodes at the current simulation timestep. See Physical quantities below for more details.

• time: a TimeInfo object containing information such as the current time and timestep;

• stats: a SimulationStats object containing information such as the total number of reconnections since the beginning of the simulation;

• to: a TimerOutput, which records the time spent on different functions;

• cache_bs: the Biot–Savart cache, which contains data from short- and long-range computations.

Physical quantities

The quantities field is a NamedTuple containing instantaneous data on filament nodes. For example, quantities.vs is the self-induced vortex velocity. To get the velocity of filament i at node j, one can do quantities.vs[i][j]. Moreover, these quantities are often interpolable, which means that one can do quantities.vs[i](j, 0.5) to get the velocity in-between two nodes.

For convenience, if one has a VortexFilamentSolver iter, the shortcut iter.vs is equivalent to iter.quantities.vs (this also applies to all other quantities).

List of quantities

• vs: self-induced vortex velocity (due to Biot–Savart law). If enabled, this also includes the contribution of an external_velocity (see init).

• ψs: self-induced streamfunction on vortices. This is useful for estimating the kinetic energy associated to the induced velocity field. If enabled, this also includes the contribution of an external_streamfunction (see init).

If a normal fluid is present (e.g. by passing forcing = NormalFluidForcing(...) to init), then the following quantities are also available:

• vL: actual filament velocity after mutual friction due to normal fluid (see NormalFluidForcing).

• vn: normal fluid velocity at vortex locations.

• tangents: unit tangents $\bm{s}'$ at vortex locations.

For convenience, if there is no normal fluid, then vL is defined as an alias of vs (they're the same object).

solve!(iter::VortexFilamentSolver)

Advance vortex filament solver to the ending time.

See also step! for advancing one step at a time.

step!(iter::VortexFilamentSolver) -> SimulationStatus

Advance solver by a single timestep.

Returns a SimulationStatus. Currently, this can be:

• SUCCESS,
• NO_VORTICES_LEFT: all vortices have been removed from the system; the simulation should be stopped.

TemporalScheme

Abstract type representing a timestepping scheme.

ExplicitScheme <: TemporalScheme

Abstract type defining an explicit temporal scheme.

Euler <: ExplicitScheme

Standard first-order Euler scheme.

Midpoint <: ExplicitScheme

Second-order, two-stage explicit midpoint method.

RK4 <: ExplicitScheme

Classic 4-stage Runge–Kutta method.

SSPRK33 <: ExplicitScheme

Three-stage, third-order strong stability preserving (SSP) method.

See Wikipedia for details.

DP5 <: ExplicitScheme

Dormand–Prince embedded 4/5 Runge–Kutta method.

For now, the implementation is experimental and can be further optimised.

See Wikipedia for details.

ImplicitExplicitScheme <: TemporalScheme

Abstract type defining an implicit-explicit (a.k.a. IMEX) temporal scheme.

The defined IMEX schemes treat the localised induction approximation (LIA) term as implicit. This may allow to increase the timestep, as it is the LIA term which imposes a small timestep in VFM simulations. Moreover, since this term is quite cheap to compute (compared to non-local interactions), treating it implicitly is not very expensive.

IMEXEuler <: ImplicitExplicitScheme

Forward-backward Euler scheme.

This is the (1,2,1) scheme in the notation of Ascher et al. (1997).

Ascher343 <: ImplicitExplicitScheme

3rd order, 4 stage IMEX Runge–Kutta scheme by Ascher et al. (Appl. Numer. Math., 1997).

⚠ This scheme may be removed in the future as it behaves very similarly to KenCarp3.

KenCarp3 <: ImplicitExplicitScheme

3rd order, 4 stage IMEX Runge–Kutta scheme by Kennedy & Carpenter (Appl. Numer. Math., 2003).

KenCarp4 <: ImplicitExplicitScheme

4th order, 6 stage IMEX Runge–Kutta scheme by Kennedy & Carpenter (Appl. Numer. Math., 2003).

SplittingScheme <: TemporalScheme

Abstract type defining a splitting scheme (such as Strang splitting).

Like IMEX (ImplicitExplicitScheme) and multirate (MultirateScheme) schemes, the idea is to split the Biot–Savart integral into two terms, respectively governing the "fast" dynamics (usually the local term) and the "slow" dynamics (non-local interactions).

The evolution due to both terms is approximated using some kind of Runge–Kutta scheme.

Strang([fast = RK4()], [slow = Midpoint()]; nsubsteps::Int = 1) <: SplittingScheme

2nd order Strang splitting scheme.

Uses one scheme for advancing the "fast" terms (assumed to be cheap to compute as well), and possibly a different scheme for the "slow" (and expensive) terms. By default these schemes are respectively taken to be RK4 and the 2nd order Midpoint method.

By default, according to Strang splitting, the fast term is advanced with a timestep of dt/2 (twice in a full timestep). One can pass nsubsteps to use even smaller timesteps for the fast term.

See SplittingScheme for more details.

MultirateScheme <: TemporalScheme

Abstract type defining a multirate scheme.

The idea is to treat different terms of the evolution equations with different timesteps (and possibly different integration schemes). Concretely, the fast dynamics is integrated with a smaller timestep (and a different inner scheme) than the slowly-evolving motions.

MultirateMidpoint([inner = Euler()], M::Int) <: MultirateScheme

2nd order, two-stage multirate infinitesimal, generalised additive Runge–Kutta (MRI-GARK) scheme.

Uses the explicit RK scheme inner for the fast component (by default a 1st order Euler scheme), with M inner steps for each outer RK stage.

This is the MRI-GARK-ERK22a method from Sandu, SIAM J. Numer. Anal. 57 (2019).

SanduMRI33a([inner = Midpoint()], M::Int) <: MultirateScheme

3rd order, 3-stage multirate infinitesimal, generalised additive Runge–Kutta (MRI-GARK) scheme.

Uses the explicit RK scheme inner for the fast component (by default a 2nd order midpoint scheme), with M inner steps for each outer RK stage.

This is the MRI-GARK-ERK33a method from Sandu, SIAM J. Numer. Anal. 57 (2019).

SanduMRI45a([inner = SSPRK33()], M::Int) <: MultirateScheme

4th order, 5-stage multirate infinitesimal, generalised additive Runge–Kutta (MRI-GARK) scheme.

Uses the explicit RK scheme inner for the fast component, with M inner steps for each outer RK stage.

This is the MRI-GARK-ERK45 method from Sandu, SIAM J. Numer. Anal. 57 (2019).

ImplicitScheme <: TemporalScheme

Abstract type representing an implicit (possibly multi-stage) scheme.

These kinds of schemes are mainly meant to be used as inner schemes (for resolving fast-evolving dynamics which are also cheap to compute) when using MultirateScheme.

CrankNicolson() <: ImplicitScheme

2nd order Crank–Nicolson implicit scheme.

LocalTerm <: FastBiotSavartTerm

Identifies fast dynamics with the local (LIA) term associated to the desingularisation of the Biot–Savart integral.

This is useful for split timestepping schemes like IMEX or multirate methods.

ShortRangeTerm <: FastBiotSavartTerm

Identifies fast dynamics with the short-range component of Ewald splitting.

This is only useful for split timestepping schemes like IMEX or multirate methods.

AdaptivityCriterion

Abstract type representing a temporal adaptivity criterion.

Implemented adaptivity criteria are:

• NoAdaptivity: disables time adaptivity;

• AdaptBasedOnSegmentLength: determines the timestep $Δt$ based on the minimum distance $ℓ_{\min}$ between two filament nodes ($Δt ∝ ℓ_{\min}^{-2}$);

• AdaptBasedOnVelocity: determines the timestep $Δt$ based on the maximum velocity $v_{\max}$ of filament nodes and on a predefined distance $δ$ ($Δt = δ / v_{\max}$).

Combining multiple criteria

Adaptivity criteria can be combined using |. Example:

adaptivity = AdaptBasedOnSegmentLength(1.4) | AdaptBasedOnVelocity(0.01)

As expected, the timestep $Δt$ will be chosen so that it satisfies both criteria.

NoAdaptivity <: AdaptivityCriterion
NoAdaptivity()

Disable temporal adaptivity, leaving the timestep $Δt$ constant.

AdaptBasedOnSegmentLength <: AdaptivityCriterion
AdaptBasedOnSegmentLength(γ::Float64)

Adapt timestep $Δt$ based on the minimum distance $δ$ between two filament nodes.

More precisely, the timestep is set to $Δt = γ \, T_{\text{KW}}(δ)$, where $γ$ is a dimensionless factor to be chosen, and:

\[T_{\text{KW}}(λ) = \frac{2 λ^2}{Γ} {\left[ \ln\left( \frac{λ}{πa} \right) + \frac{1}{2} - (Δ + γ) \right]}^{-1}\]

is the period of a Kelvin wave of wavelength $λ$. See ParamsBiotSavart for the definitions of the vortex parameters $Γ$, $a$ and $Δ$.

This criterion is somewhat analogous to the CFL condition in grid-based computations, and $γ$ is the analogous of the maximum CFL number to be allowed. As such, the right value of $γ$ will depend on the chosen temporal scheme.

For example, the RK4 scheme seems to require $γ ≈ 1$ to remain stable.

AdaptBasedOnVelocity <: AdaptivityCriterion
AdaptBasedOnVelocity(δ_crit::Float64; safety_factor = 0.8)

Adapt timestep $Δt$ based on the maximum velocity $v_{\max}$ of filament nodes.

The objective is that the resulting maximum displacement $δ_{\max} = v_{\max} Δt$ stays below the given critical displacement δ_crit.

One application of this criterion is to ensure that reconnections happen in-between two solver iterations (that is, to avoid that two filaments cross each other without reconnecting). In this case, $δ$ should be proportional to the chosen distance below which reconnections are performed.

Implementation details

This criterion is used in two ways:

1. A priori, to decide the $Δt$ to be used in the next iteration given the velocities of filament nodes at the start of the iteration (at time $t$). This ensures $Δt ≤ δ_{\text{crit}} / v_{\max}$ where $v_{\max}$ is computed at the start of the iteration.

2. A posteriori, after the actual node displacements $δ$ from time $t$ to time $t + Δt$ have been computed (e.g. using some Runge–Kutta scheme). If some $δ$ is larger than $δ_{\text{crit}}$, then the iteration is recomputed after halving the timestep ($Δt → Δt/2$). In this case, the original displacements are thrown away (rejected). This process can be eventually repeated until the criterion is satisfied.

Optional safety factor

The optional safety_factor should be ≤1, which will further reduce the timestep chosen a priori in step 1. This is to try to avoid too many rejected timesteps in step 2, e.g. in case the actual advection velocity is larger than the velocity at the beginning of the timestep.

In combination with other criteria

In principle, using this criterion can lead to an infinite timestep when the velocities are zero. For this reason, it's a good idea to combine this criterion with the AdaptBasedOnSegmentLength or the MaximumTimestep criterion. For example:

adaptivity = AdaptBasedOnVelocity(2.0) | AdaptBasedOnSegmentLength(0.9)
adaptivity = AdaptBasedOnVelocity(2.0) | MaximumTimestep(0.01)

In fact, the second option is done automatically in init if only an AdaptBasedOnVelocity is passed. In that case, the maximum timestep is taken to be the dt passed to init.

MaximumTimestep <: AdaptivityCriterion
MaximumTimestep(Δt_max::Float64)

Criterion ensuring that the timestep will be kept below a maximal value Δt_max.

inject_filament!(iter::VortexFilamentSolver, f::AbstractFilament)

Inject a filament onto a simulation.

This function must be called from within the affect! function attached to the solver (see init for details). Otherwise, the solver will likely use garbage values for the velocity of the injected filament, leading to wrong results or worse.

Basic usage

A very simplified (and incomplete) example of how to inject filaments during a simulation:

# Define function that will be called after each simulation timestep
function my_affect!(iter::VortexFilamentSolver)
    f = Filaments.init(...)  # create new filament
    inject_filament!(iter, f)
    return nothing
end

# Initialise and run simulation
prob = VortexFilamentProblem(...)
iter = init(prob, RK4(); affect! = my_affect!, etc...)
solve!(iter)

TemporalSchemeCache{Scheme <: TemporalScheme}

Contains buffers needed by a temporal scheme.

TimeInfo

Contains information on the current time and timestep of a solver.

Some useful fields are:

• t::Float64: current time;

• dt::Float64: timestep to be used in next iteration;

• dt_prev::Float64 : timestep used in the last performed iteration;

• nstep::Int: number of timesteps performed until now;

• nrejected::Int: number of rejected iterations.

When using the AdaptBasedOnVelocity criterion, an iteration can be rejected if the actual filament displacement is too large compared to what is imposed by the criterion. In that case, the iteration will be recomputed with a smaller timestep (dt → dt/2).

SimulationStats{T <: AbstractFloat}

Contains accumulated statistics of different events occurring since the beginning of a simulation.

Available fields

• reconnection_count::Int: total number of reconnections;

• reconnection_length_loss::T: accumulated decrease of filament length due to reconnections. This is estimated from each reconnection event as the difference between the local vortex lengths before and after the reconnection, using the straight segment approximation.

• filaments_removed_count::Int: total number of removed filaments (this generally happens when the number of discretisation points becomes too small for the spatial discretisation to work);

• filaments_removed_length::T: total length of removed filaments. Note that this length is estimated using a straight segment approximation (no quadratures). This is because filaments are removed when they can no longer be represented using a continuous interpolation function.

maximum_displacement(::AdaptivityCriterion) -> Float64

Return the maximum node displacement $δ_{\text{crit}}$ allowed by the adaptivity criterion.

More precisely, $δ_{\text{crit}}$ controls the maximum displacement of a filament node during a single timestep of duration $Δt$.

This function is mainly relevant when using the AdaptBasedOnVelocity criterion. For other criteria this returns Inf, meaning that they don't limit the maximum allowed displacement.