API Reference

Types

Kernels

ContourDynamics.EulerKernel Type

EulerKernel <: AbstractKernel

Kernel for 2-D Euler (barotropic) vortex dynamics, using the Green's function G(r) = -log(r) / (2π).

source

ContourDynamics.QGKernel Type

QGKernel{T}(Ld)

Kernel for single-layer quasi-geostrophic dynamics with deformation radius Ld. Uses the modified Bessel function K₀(r/Ld) as the Green's function.

source

ContourDynamics.SQGKernel Type

SQGKernel{T}(delta)

Kernel for surface quasi-geostrophic (SQG) dynamics using the fractional Laplacian Green's function G(r) = -1/(2πr).

The SQG velocity at a contour boundary is singular (the tangential component diverges logarithmically), so a regularization length delta > 0 is required. The kernel 1/r is replaced by 1/√(r² + δ²). A typical choice is delta ≈ μ, the minimum segment length used for surgery.

Sign convention

Positive PV (surface buoyancy) induces counter-clockwise circulation, matching the EulerKernel convention.

source

ContourDynamics.MultiLayerQGKernel Type

MultiLayerQGKernel{N,M,T}(Ld, coupling)

Kernel for N-layer quasi-geostrophic dynamics with M = N-1 deformation radii Ld and layer coupling matrix coupling. The constructor eigen-decomposes the coupling matrix for efficient velocity evaluation.

The coupling matrix should already incorporate layer thicknesses (i.e. it is the full stretching operator, not raw interface stretching coefficients).

Note

The evolution and energy diagnostics derive modal deformation radii from the coupling matrix eigenvalues (Ld_mode = 1/√|λ|). The constructor therefore validates that the provided Ld values match the nonzero modal radii implied by coupling, preventing silent simulation of a different physical model.

source

Contours and Domains

ContourDynamics.PVContour Type

PVContour{T}(nodes, pv[, wrap])

A piecewise-linear contour carrying a potential-vorticity jump pv. nodes is the ordered list of vertices. For contours that span a periodic domain, wrap gives the shift vector that closes the final segment back to nodes[1]; it defaults to zero for ordinary closed contours.

source

ContourDynamics.UnboundedDomain Type

UnboundedDomain <: AbstractDomain

An infinite, unbounded two-dimensional domain.

source

ContourDynamics.PeriodicDomain Type

PeriodicDomain{T}(Lx, Ly)

A doubly-periodic rectangular domain with half-widths Lx and Ly, i.e. the domain [-Lx, Lx) × [-Ly, Ly).

source

Problem Structs

ContourDynamics.ContourProblem Type

ContourProblem{K,D,T,Dev}(kernel, domain, contours; dev=CPU())

A single-layer contour-dynamics problem with a velocity kernel, computational domain, a vector of PVContours, and a target device (CPU or GPU).

source

ContourDynamics.MultiLayerContourProblem Type

MultiLayerContourProblem{N,K,D,T,Dev}(kernel, domain, layers; dev=CPU())

An N-layer contour-dynamics problem. Each element of the layers tuple holds the contours for one layer. The optional dev keyword selects the target device (CPU or GPU) for buffer allocation.

source

ContourDynamics.SurgeryParams Type

SurgeryParams{T}(delta, mu, Delta_max, area_min, n_surgery)

Parameters controlling contour surgery.

• delta: proximity threshold for detecting close contour segments.

• mu: minimum segment length after remeshing.

• Delta_max: maximum segment length after remeshing.

• area_min: minimum enclosed area; contours smaller than this are removed.

• n_surgery: number of time-steps between surgery passes.

Leapfrog stepper interaction

Surgery may change the number of nodes (via reconnection, filament removal, and remeshing), which invalidates the previous-step history required by LeapfrogStepper. After each surgery pass the leapfrog method re-bootstraps with an RK2 midpoint half-step, maintaining second-order accuracy. If high-order accuracy is important, prefer RK4Stepper or increase n_surgery to reduce the frequency of re-bootstrapping.

source

ContourDynamics.Problem Type

Problem{P,S,SP}

Convenience wrapper bundling a ContourProblem (or MultiLayerContourProblem), a time stepper, and optional SurgeryParams into a single object.

Construct via the keyword factory Problem(; kwargs...) or directly:

Problem(contour_problem, stepper, surgery_params)

source

Accessors

ContourDynamics.contours Function

Return the contours of the underlying problem.

source

ContourDynamics.kernel Function

Return the kernel of the underlying problem.

source

ContourDynamics.domain Function

Return the domain of the underlying problem.

source

Time Steppers

ContourDynamics.RK4Stepper Type

RK4Stepper{T,A}(dt, n; dev=CPU())

Classical fourth-order Runge–Kutta time stepper with step size dt. Allocates internal buffers for n nodes on the given device.

source

ContourDynamics.LeapfrogStepper Type

LeapfrogStepper{T,A}(dt, n; dev=CPU(), ra_coeff=0.05)

Leapfrog (second-order centred) time stepper with step size dt. The first step is bootstrapped with an RK2 midpoint half-step. Buffers are allocated on the given device.

source

Velocity Computation

ContourDynamics.velocity! Function

velocity!(vel, prob::ContourProblem)

Compute velocity at every contour node of prob, storing results in vel. Uses the proxy-FMM path only when that experimental accelerator is explicitly enabled. Otherwise, large single-layer problems use the production treecode path and small problems fall back to the validated direct evaluator.

source

velocity!(vel, prob::MultiLayerContourProblem)

Compute velocity at all nodes across all layers using modal decomposition. Uses the treecode for large problems, the proxy FMM when explicitly enabled, and direct O(N²) evaluation for small problems.

source

ContourDynamics.velocity Function

velocity(prob::ContourProblem, x::SVector{2,T})

Compute velocity at a single point x from all contours in prob.

source

ContourDynamics.segment_velocity Function

segment_velocity(::QGKernel, ::UnboundedDomain, x, a, b)

Velocity at point x due to a vortex patch contour segment from a to b using the QG Green's function G(r) = K₀(r/Ld) / (2π).

The contour dynamics velocity is: v_seg = (1/(2π)) ∫₀¹ K₀(|x-P(t)|/Ld) ds dt

Uses singular subtraction: the log singularity in K₀ is handled analytically (matching the Euler kernel), and the smooth remainder [K₀(r/Ld) + log(r)] is integrated with 5-point Gauss-Legendre quadrature.

source

segment_velocity(::SQGKernel, ::UnboundedDomain, x, a, b)

Velocity at point x due to a surface buoyancy patch contour segment from node a to node b with unit buoyancy jump, using the regularized SQG Green's function G(r) = -1/(2π√(r²+δ²)).

The contour integral is: v_seg = -(1/(2π)) t̂ [F(u_a) - F(u_b)]

where F(u) = log(u + √(u² + h_eff²)) and h_eff² = h² + δ².

The -(1/(2π)) prefactor (vs -(1/(4π)) for Euler) reflects the different Green's function normalisation: SQG uses 1/(2πr) while Euler uses log(r²)/(4π).

source

segment_velocity(kernel::EulerKernel, domain::PeriodicDomain, x, a, b)

Velocity at point x from segment a→b in a periodic domain using Ewald summation.

Uses singular subtraction: the log singularity is handled analytically by the unbounded Euler segment velocity, and only the smooth periodic correction G_per - G_∞ is integrated with 3-point Gauss-Legendre quadrature.

The periodic correction decomposes as:

• Central-image real-space: (1/(4π))[E₁(α²r²) + log(r²)] → (1/(4π))(-γ - 2ln α) as r→0

• Non-central real-space: (1/(4π)) Σ_{images≠0} E₁(α²|r+shift|²) (smooth)

• Fourier space: Σ_{k≠0} coeff * cos(k·r) (smooth)

source

segment_velocity(kernel::QGKernel, domain::PeriodicDomain, x, a, b)

Velocity at point x from segment a→b in a periodic domain using the QG kernel.

Decomposes the periodic QG Green's function as: G_QG_per = G_Euler_per - G_correction where the Euler part is handled by the validated Ewald summation, and the correction is a smooth, rapidly convergent Fourier series: G_corr(r) = -(1/A) Σ_{k≠0} cos(k·r) κ²/(k²(k²+κ²)) with κ = 1/Ld. Coefficients decay as 1/k⁴, so the truncated sum converges without Gaussian damping.

source

segment_velocity(kernel::SQGKernel, domain::PeriodicDomain, x, a, b)

Velocity at point x from segment a→b in a periodic domain using the SQG kernel.

Uses singular subtraction: the regularized unbounded SQG velocity (exact arcsinh antiderivative with $\delta >0$) handles the $1/r$ singularity analytically, and the smooth periodic correction $G_{\text{per}}-G_{\mathrm{\infty }}$ is integrated with 3-point Gauss-Legendre quadrature.

The periodic correction decomposes as:

• Central-image real-space: $-(1/(2\pi ))[\mathrm{erfc}(\alpha r)/r-1/\sqrt{r^2+{\delta }^2}]$ (finite at GL quadrature points since $r>0$)

• Non-central real-space: $-(1/(2\pi ))\mathrm{erfc}(\alpha |\mathbf{r}-{\mathbf{R}}_n|)/|\mathbf{r}-{\mathbf{R}}_n|$

• Fourier space: $-(1/(2\pi ))\sum c_k\cos (\mathbf{k}\cdot \mathbf{r})$ with $c_k=(2\pi /|k|)e^{-k^2/(4{\alpha }^2)}/A$

source

Time Integration

ContourDynamics.timestep! Function

timestep!(prob::ContourProblem, stepper::RK4Stepper)

Advance all contour nodes by one RK4 step.

source

timestep!(prob::ContourProblem, stepper::LeapfrogStepper)

Advance all contour nodes by one leapfrog step. First step uses RK2 midpoint method to bootstrap.

source

ContourDynamics.evolve! Function

evolve!(prob::Problem; nsteps, callbacks=nothing)

Run the simulation for nsteps time steps. Surgery is applied according to prob.surgery_params (or skipped if nothing).

source

ContourDynamics.resize_buffers! Function

resize_buffers!(stepper::RK4Stepper, prob::ContourProblem)

Resize RK4 work arrays after surgery changes node count. Stepper buffers are always CPU Vector (even for GPU problems), since the GPU velocity path handles its own device allocation internally.

source

resize_buffers!(stepper::LeapfrogStepper, prob::ContourProblem)

Resize leapfrog work arrays after surgery. Resets initialization flag since node correspondence is lost. Stepper buffers are always CPU Vector (even for GPU problems).

source

Surgery

ContourDynamics.surgery! Function

surgery!(prob::ContourProblem, params::SurgeryParams)

Full Dritschel surgery suite: remesh → reconnect → remove filaments. Mutates prob.contours in place.

source

ContourDynamics.remesh Function

remesh(c::PVContour, params::SurgeryParams; _buf=nothing, _arc_buf=nothing, _vnodes_buf=nothing)

Redistribute nodes along contour c so that every segment length lies between params.mu and params.Delta_max. Returns a new PVContour.

The optional _buf, _arc_buf, and _vnodes_buf keywords accept pre-allocated vectors that are reused across calls to avoid repeated heap allocation (internal optimisation used by surgery!).

source

ContourDynamics.find_close_segments Function

find_close_segments(contours, spatial_index, delta[, domain])

Find pairs of contour segments whose closest approach is within delta, using the spatial index for candidate filtering. Returns vector of (ci, i, cj, j) tuples where i,j are segment indices (each segment goes from node i to next_node(c, i)).

For PeriodicDomain, minimum-image distances are used so that segments close across the periodic boundary are correctly detected.

source

ContourDynamics.build_spatial_index Function

Build a spatial index for all contour segments, binned by grid of size delta. Each segment is binned at its endpoint and at its midpoint so that long segments whose interiors cross a bin boundary are discoverable via neighbour-bin queries.

For PeriodicDomain, coordinates are wrapped to the canonical domain and ghost entries are inserted near boundaries so that segments physically close across the periodic seam appear in adjacent bins.

source

ContourDynamics.reconnect! Function

reconnect!(contours, close_pairs[, domain])

Perform contour reconnection for identified close segment pairs. Same contour → split; different contours with same PV → merge.

Each sub-contour produced by a split contains both pinch-point nodes as its first and last vertices, ensuring a well-formed closing segment. Merged contours are stitched so that traversal orientation is preserved.

Warning

Reconnection produces near-duplicate nodes at stitch points. Callers should remesh all contours after reconnection to clean up short/long segments. The top-level surgery! does this automatically.

source

ContourDynamics.remove_filaments! Function

Remove contours with |area| < area_min. Spanning contours are always kept.

source

Diagnostics

ContourDynamics.vortex_area Function

vortex_area(c::PVContour)

Signed area enclosed by contour c using the shoelace formula.

source

ContourDynamics.centroid Function

centroid(c::PVContour)

Centroid of the region enclosed by contour c, via Green's theorem.

source

ContourDynamics.ellipse_moments Function

ellipse_moments(c::PVContour)

Second moments → (aspect_ratio, orientation_angle).

source

ContourDynamics.energy Function

energy(prob)

Kinetic energy of a ContourProblem or MultiLayerContourProblem, computed via contour integrals of the appropriate Green's function.

source

ContourDynamics.enstrophy Function

enstrophy(prob)

Enstrophy ½ ∑ qᵢ² Aᵢ of a ContourProblem or MultiLayerContourProblem.

Uses signed area Aᵢ from the shoelace formula (positive for CCW, negative for CW contours), so contributions from inner boundaries are subtracted. Spanning contours are excluded (see circulation).

Warning

This diagnostic is exact when contours encode disjoint PV regions through their signed areas, but it does not reconstruct the fully squared piecewise PV field for arbitrary nested multi-jump contour sets. In those cases the missing cross-terms make the result only approximate.

source

ContourDynamics.circulation Function

circulation(prob)

Total circulation Γ = ∑ qᵢ Aᵢ of a ContourProblem or MultiLayerContourProblem.

Warning

Spanning contours have undefined area and are silently excluded. If your problem uses spanning contours (e.g. from beta_staircase), the returned circulation only reflects closed contours.

source

ContourDynamics.angular_momentum Function

angular_momentum(prob)

Angular momentum ∑ qᵢ ∫ r² dA of a ContourProblem or MultiLayerContourProblem.

source

Contour Helpers

ContourDynamics.nnodes Function

nnodes(c::PVContour)

Number of nodes (vertices) in contour c.

source

ContourDynamics.nlayers Function

nlayers(kernel_or_problem)

Return the number of layers in a MultiLayerQGKernel or MultiLayerContourProblem.

source

ContourDynamics.total_nodes Function

total_nodes(prob)

Total number of nodes across all contours in a ContourProblem or MultiLayerContourProblem.

source

ContourDynamics.arc_lengths Function

arc_lengths(c::PVContour)

Return a vector of segment lengths for each consecutive node pair in contour c.

source

ContourDynamics.next_node Function

Get the next node after index j, handling periodic wrap for spanning contours.

source

ContourDynamics.beta_staircase Function

beta_staircase(beta, domain::PeriodicDomain, n_steps; T=Float64)

Create spanning contours that discretize the background PV gradient βy into a PV staircase on a periodic domain.

Returns a Vector{PVContour{T}} of horizontal spanning contours, each carrying PV jump Δq = β * Δy where Δy = 2*Ly / n_steps. Nodes run left-to-right at evenly spaced y-levels from -Ly + Δy to Ly - Δy (excluding boundaries).

The wrap field is set to (2*Lx, 0) so the closing segment connects the rightmost node back to the leftmost node shifted by one period.

source

ContourDynamics.is_spanning Function

is_spanning(c::PVContour) -> Bool

Return true if contour c spans the periodic domain (i.e. has a non-zero wrap vector).

source

Shape Helpers

ContourDynamics.circular_patch Function

circular_patch(R, N, pv; cx=0.0, cy=0.0, T=Float64)

Create a circular PVContour with radius R, N nodes, and PV jump pv, centered at (cx, cy).

source

ContourDynamics.elliptical_patch Function

elliptical_patch(a, b, N, pv; cx=0.0, cy=0.0, θ=0.0, T=Float64)

Create an elliptical PVContour with semi-axes a and b, N nodes, PV jump pv, centered at (cx, cy), rotated by angle θ (radians).

source

ContourDynamics.rankine_vortex Function

rankine_vortex(R, N, Γ; cx=0.0, cy=0.0, T=Float64)

Create a Rankine vortex (uniform vorticity patch) with radius R, N nodes, and total circulation Γ, centered at (cx, cy). Returns a single-element Vector{PVContour{T}} with pv = Γ / (π R²).

source

Ewald Summation

ContourDynamics.EwaldCache Type

EwaldCache{T}

Precomputed data for Ewald summation in periodic domains.

source

ContourDynamics.build_ewald_cache Function

build_ewald_cache(domain::PeriodicDomain, kernel::EulerKernel; n_fourier=8, n_images=2)

Precompute Fourier-space coefficients for Ewald summation.

source

build_ewald_cache(domain::PeriodicDomain, kernel::SQGKernel; n_fourier=8, n_images=2)

Ewald cache for SQG kernel in periodic domain.

The SQG Green's function G(r) = -1/(2πr) is split via Ewald summation. Fourier coefficients are (2π/|k|) exp(-k²/(4α²)) / A, reflecting the fractional Laplacian's half-order (1/|k| vs Euler's 1/k²).

source

build_ewald_cache(domain::PeriodicDomain, kernel::QGKernel; n_fourier=8, n_images=2)

Ewald cache for QG kernel in periodic domain.

source

ContourDynamics.setup_ewald_cache! Function

setup_ewald_cache!(domain, kernel; n_fourier=8, n_images=2)

Pre-build and store an Ewald cache with custom parameters. Call this before evolve! to override the default n_fourier=8, n_images=2. The cached result is used automatically by all subsequent velocity computations on the same domain/kernel combination.

source

setup_ewald_cache!(domain, kernel::QGKernel; n_fourier=8, n_images=2)

Pre-build Ewald caches for QG velocity computation. The periodic QG velocity decomposes as G_QG_per = G_Euler_per - G_correction, so the velocity path uses an Euler Ewald cache internally. This method builds both the QG-specific cache and the Euler cache that the velocity path actually reads, ensuring that custom n_fourier/n_images parameters take effect.

source

ContourDynamics.clear_ewald_cache! Function

Clear all cached Ewald data.

source

Periodic Domains

ContourDynamics.wrap_nodes! Function

wrap_nodes!(prob::ContourProblem{K, PeriodicDomain{T}})

Wrap all non-spanning contour nodes into the fundamental domain. Spanning contours are left untouched since their positions encode the cross-domain topology via the wrap vector.

source

No-op for unbounded domains — nodes don't need wrapping.

source

Device Types

ContourDynamics.CPU Type

CPU()

Select CPU execution. This is the default device.

source

ContourDynamics.GPU Type

GPU()

Select GPU execution. Requires using CUDA to activate the CUDA extension.

source

ContourDynamics.device_array Function

device_array(dev)

Return the array constructor for the given device: Array for CPU(), CuArray for GPU() (when CUDA extension is loaded).

source

Internals

ContourDynamics._collect_all_nodes! Function

_collect_all_nodes!(buf, prob::ContourProblem)

Collect all nodes into pre-allocated buffer buf (in-place, non-allocating).

source

ContourDynamics._scatter_nodes! Function

_scatter_nodes!(prob::ContourProblem, all_nodes)

Write flat node vector back into contour node arrays.

source

ContourDynamics._scatter_shifted! Function

_scatter_shifted!(prob, base, delta, scale)

Write base[i] + scale * delta[i] into contour nodes without allocating.

source

ContourDynamics._expint_e1 Function

_expint_e1(x)

Compute the exponential integral E₁(x) = ∫_x^∞ e^{-t}/t dt for x > 0.

source

ContourDynamics._segment_min_dist2 Function

_segment_min_dist2(a1, b1, a2, b2)

Minimum squared distance between segments a1→b1 and a2→b2.

Uses the full parametric closest-point algorithm (Ericson, "Real-Time Collision Detection") to handle all cases: endpoint-to-segment, endpoint-to-endpoint, and interior-to-interior closest points.

source

ContourDynamics._best_stitch_nodes Function

_best_stitch_nodes(c1, i1, c2, i2[, domain])

Given that segment i1 of contour c1 is close to segment i2 of contour c2, find the pair of node indices (one from each contour) that are closest. Returns (best_i1, best_i2) — node indices into c1.nodes and c2.nodes.

For PeriodicDomain, minimum-image distances are used.

source

ContourDynamics._check_spanning_proximity Function

_check_spanning_proximity(contours, delta[, domain])

Warn if any closed contour node is within delta of a spanning contour node. This situation cannot be resolved by surgery (spanning contours are exempt from reconnection) and may indicate insufficient resolution or an overly large delta.

For PeriodicDomain, minimum-image distances are used.

source