Operators

nonlinear!(nonlin, u, ω)

Compute the nonlinear advection term in-place.

This function updates the nonlinear term array nonlin by evaluating the nonlinear contribution at every grid point, based on the velocity field u and the vorticity field ω. It represents the convective term the convective term (u · ∇)u = u × ω for incompressible flow.

The computation is parallelized over the grid using the appropriate backend (e.g. CPU or GPU).

Arguments

• nonlin: array (or array of arrays) storing the nonlinear term; modified in place.
• u: velocity field.
• ω: vorticity field.

Returns

The updated nonlin field.

nonlinear(i, u, ω, I)

Compute the nonlinear advection term for component i at grid point I.

This function evaluates the local contribution of the nonlinear term the i-th component of the cross product u × ω computed using bilinear interpolation of the velocity and vorticity fields. velocity and vorticity fields. It is called internally by nonlinear!.

Arguments

• i: index of the velocity component being computed.
• u: velocity field.
• ω: vorticity field.
• I: Cartesian grid index.

Returns

The scalar nonlinear term at the specified component and grid location.

rot!(ω, u; h)

Compute the vorticity field ω from the velocity field u in-place.

Arguments

• ω: Array of arrays where the computed vorticity components will be stored (mutated in-place).
• u: Array of arrays representing the velocity field.
• h: Grid spacing used for finite-difference approximation of the curl (keyword argument).

Description

This function updates ω directly by looping over each component and grid index, calling rot(i, u, I; h) for each point. The computation is parallelized using the available backend (CPU or GPU) for efficiency.

Returns

• ω: The updated vorticity field (same array as input, modified in-place).

rot(i, u, I; h)

Compute the i-th component of the vorticity (curl) at a single grid point I from a velocity field u.

Arguments

• i: Index of the vorticity component to compute (e.g., 1 for x, 2 for y).
• u: Array of arrays representing the velocity field.
• I: Cartesian index of the grid point where the curl is computed.
• h: Grid spacing used for finite-difference approximation.

Returns

• Scalar value representing the i-th component of the vorticity at point I.

Notes

• This function computes the curl at a single point. To compute the full vorticity field over the grid, use rot!(ω, u; h), which applies rot in-place across all grid points.
• The function uses finite differences and cross-product indexing (sumcross) to compute the curl:

(∇ × u)i = sum{(j,k)} (uk[I] - uk[I-δ(j)]) / h.

curl!(u, ψ; h)

Compute the velocity field u as the curl of a potential field ψ over the entire grid.

Arguments

• u: Array of arrays representing the velocity field (updated in-place).
• ψ: Array of arrays representing the potential field.
• h: Grid spacing used for finite-difference approximation.

Returns

• The updated velocity field u.

Notes

• This function applies the curl computation at every grid point in-place.
• It uses finite differences and the helper function curl(i, ψ, I; h) to compute each component.
• The backend is automatically chosen for parallel or GPU execution.

curl(i, ψ, I; h)

Compute the i-th component of a velocity field as the curl of a scalar potential ψ at a specific grid point.

Arguments

• i: Index of the velocity component to compute (e.g., 1 for x, 2 for y).
• ψ: Array of arrays representing the scalar potential field.
• I: Cartesian index of the grid point where the curl is evaluated.
• h: Grid spacing used for finite-difference approximation.

Returns

• A scalar value representing the i-th component of the curl at grid point I.

Notes

• The result is divergence-free by construction.
• Uses centered finite differences and sumcross to handle component indices.

struct LaplacianPlan

Holds the data required to compute the Laplacian efficiently in spectral space using FFTs.

Fields

• λ: eigenvalues of the Laplacian for spectral multiplication.
• work: temporary workspace array.
• fwd: forward FFT plan.
• inv: inverse FFT plan.
• n_logical: logical size of the FFT domain.

Constructor

LaplacianPlan(ω_i, i, n)

Creates a LaplacianPlan for a given field component ω_i on a grid of size n (e.g., SVector(nx, ny)).

Arguments

• ω_i: array representing one component of the field (e.g., vorticity).
• i: component index (1, 2, or 3 for x, y, z).
• n: grid resolution as SVector{N}.

Returns

• A LaplacianPlan struct containing all data needed to:
  1. Transform the field to spectral space.
  2. Multiply by Laplacian eigenvalues.
  3. Transform back to physical space.

Notes

• The constructor sets up FFT plans (fwd and inv) and eigenvalues (λ) automatically.
• The Laplacian is then applied as Δu ≈ λ .* fwd(u).

laplacian_fft_kind(i, nd)

Return a tuple specifying the FFT type to use along each dimension of a multidimensional array for the Laplacian operator.

Arguments

• i: The dimension index that should use a cosine transform (FFTW.REDFT01).
• nd: Total number of dimensions of the array.

Returns

• A tuple of length nd with FFTW.REDFT01 for dimension i and FFTW.RODFT00 for all other dimensions.

Notes

• REDFT01 (DCT-I) is typically used for Neumann or periodic boundary conditions.
• RODFT00 (DST-I) is typically used for Dirichlet (zero-value) boundary conditions.

laplacian_eigvals!(λ, i)

Compute the eigenvalues of the discrete Laplacian on a multidimensional grid and store them in λ.

Arguments

• λ: Array to hold the eigenvalues (modified in-place)
• i: The dimension index along which the cosine transform (DCT) is applied; other dimensions use sine transforms (DST)

Returns

• The array λ filled with Laplacian eigenvalues

Notes

• The computation runs in parallel on the backend associated with λ
• Eigenvalues correspond to the discrete Laplacian under DCT (dimension i) and DST (other dimensions)

laplacian_plans(ω, n)

Create a Laplacian plan for each component of a vector field ω.

Arguments

• ω: Vector field whose components require Laplacian plans
• n: Grid resolution (e.g., tuple of (nx, ny) or (nx, ny, nz))

Returns

• A tuple of LaplacianPlan objects, one for each component of ω

Notes

• Each LaplacianPlan contains the eigenvalues and FFT plans needed for spectral Poisson solvers
• Useful when different components require different FFT types due to mixed boundary conditions

EigenbasisTransform

Represents a spectral transform that applies a function f (e.g., inverse Laplacian) in the eigenbasis of the Laplacian. It uses precomputed Laplacian plans to efficiently apply forward and inverse FFTs, multiply by f(λ) in spectral space, and transform back.

Fields

• f: Function or scalar applied to the eigenvalues.
• plan: Tuple of LaplacianPlan objects wrapped in an OffsetTuple for staggered grid alignment.

Constructor

• EigenbasisTransform(f, plans::Tuple) wraps a regular tuple of LaplacianPlan objects into an OffsetTuple automatically. This allows you to pass a standard tuple without needing to manually construct an OffsetTuple.

Call Overloads

• (X::EigenbasisTransform)(y, ω) applies the transform to all components of ω, storing the result in y.
• (X::EigenbasisTransform)(yᵢ, ωᵢ, i) applies the transform to the i-th component using its corresponding Laplacian plan.

Notes

• Useful for spectral Poisson solvers and other operations in the Laplacian eigenbasis.
• The OffsetTuple ensures proper alignment for staggered grid variables.

_coarse_indices(n::NTuple{N}, loc::Edge{Dual}) where {N}

Compute the index ranges for a coarser grid that is a subset of the original grid.

Description

This function determines which indices to retain when coarsening a field defined on edges of the dual grid. It is used in multigrid methods to identify the portion of the grid that remains after reducing resolution.

Arguments

• n: Tuple giving the number of grid points in each dimension.
• loc: Specifies the edge location in the dual grid.

Returns

A tuple of index ranges corresponding to the coarse grid coordinates.

multidomain_coarsen!(ω², ω¹; n)

Coarsen a fine-grid field ω¹ into a coarser field ω² for use in multigrid solvers.

Description

This function reduces the resolution of each component of a vector field by averaging fine-grid values into the coarser grid using a restriction stencil. It runs in parallel across the computational backend and relies on _coarse_indices and multidomain_coarsen for index mapping and averaging.

Arguments

• ω²: Output field on the coarser grid.
• ω¹: Input field on the finer grid.
• n: Tuple giving the grid resolution.

Returns

The coarsened field ω².

multidomain_coarsen(i, ωᵢ, I²; n)

Compute the coarse-grid value of a field component from its fine-grid representation.

Description

For component i, this function calculates the value at coarse-grid index I² by averaging corresponding fine-grid points using a predefined stencil. It is used internally by multidomain_coarsen! to define the restriction operation.

Arguments

• i: Component index.
• ωᵢ: Fine-grid array for the given component.
• I²: Cartesian index on the coarse grid.
• n: Tuple giving the grid resolution.

Returns

The coarse-grid value at index I².

_coarsen_stencil(T)

Return a normalized 3×3 smoothing stencil used for grid coarsening.

Description

This function defines a weighted averaging kernel for use in multigrid coarsening or smoothing operations. It produces a static 3×3 matrix from StaticArrays.jl with higher weight at the center and symmetric lower weights around it. The entries sum to 1, ensuring mass conservation during coarsening.

Arguments

• T: Element type of the output matrix.

Returns

A 3×3 static matrix of type T used as a normalized coarsening stencil.

_fine_indices(i, n, I)

Return the fine-grid indices corresponding to a coarse-grid index I.

Description

This function identifies which fine-grid cells map to a given coarse-grid location, used for coarsening or restriction in multigrid solvers. It supports both 2D and 3D cases:

• In 2D, it returns a single CartesianIndices range.
• In 3D, it returns two CartesianIndices planes associated with the selected component i.

Arguments

• i: Component index (used only in 3D).
• n: Tuple representing the grid size.
• I: Coarse-grid index tuple.

Returns

A tuple of one or more CartesianIndices objects representing fine-grid regions corresponding to the coarse-grid index I.

_fine_range(n, I)

Return the range of fine-grid indices corresponding to a coarse-grid index I along one dimension.

Description

Maps a 1D coarse-grid index to the corresponding three-point region on the fine grid, assuming a 4:1 refinement ratio. The range is centered around 2 * (I - n/4) and includes the neighboring points at offsets -1, 0, and +1.

Arguments

• n: Grid size along the considered dimension.
• I: Coarse-grid index.

Returns

A UnitRange of three fine-grid indices corresponding to the coarse-grid index I.

multidomain_interpolate!(ωb, ω; n)

Interpolate a fine-grid field ω onto boundary or ghost regions, storing results in ωb.

Description

This function computes interpolated boundary values from a fine-grid field ω and writes them into the corresponding boundary arrays ωb. It supports both 2D and 3D cases and calls the helper function [multidomain_interpolate] for the actual interpolation logic.

Arguments

• ωb: Output array for interpolated boundary values.
• ω: Input fine-grid field (e.g., vorticity).
• n: Grid size tuple, used for mapping fine-to-coarse indices.

Returns

The updated ωb containing interpolated boundary values.

multidomain_interpolate(ωᵢ, (i, j, k), dir, I¹; n)

Interpolate values from a fine grid onto a coarser grid or boundary face.

Description

Defines the interpolation rule used by [multidomain_interpolate!]. Two methods are available:

• 2D version (CartesianIndex{2}): performs simple linear interpolation along one direction.
• 3D version (CartesianIndex{3}): performs bilinear interpolation on a 2D plane, adjusting for component offsets.

Arguments

• ωᵢ: Component i of the fine-grid field.
• (i, j, k): Index tuple defining component and orientation.
• dir: Direction index for boundary interpolation.
• I¹: Fine-grid index at which interpolation is performed.
• n: Grid size tuple.

Returns

Interpolated scalar value for the target coarse-grid or boundary location.

set_boundary!(ω, ωb)

Copy boundary values from a boundary buffer ωb into the main field ω.

Description

This function updates the field ω by setting its boundary regions to the corresponding values stored in ωb. It is typically used after interpolation or restriction steps in multigrid or multidomain solvers to enforce consistent boundary conditions.

The operation is performed in parallel using backend-agnostic loops.

Arguments

• ω: Main field to be updated (e.g., vorticity or velocity component arrays).
• ωb: Boundary buffer containing precomputed boundary values for each component of ω.

Returns

The modified field ω with updated boundary values.

add_laplacian_bc!(Lψ, factor, ψb)

Apply boundary condition corrections to the Laplacian of a vector field.

Description

This function modifies the Laplacian Lψ of a field by adding contributions from boundary values in ψb. It ensures that Dirichlet or Neumann boundary conditions are correctly enforced in the discrete Laplacian.

Arguments

• Lψ: Vector of arrays representing the Laplacian of the field ψ.
• factor: Scalar multiplier for the boundary corrections.
• ψb: Boundary buffer holding Dirichlet or Neumann values at the domain boundaries.

Returns

The modified Lψ with boundary corrections applied (in-place).

laplacian_bc_ii(ψb, i, dir, I)

Compute the diagonal Laplacian boundary correction for component i.

Description

This function calculates the correction term for the Laplacian at a boundary along direction dir for the i-th component. It uses differences of precomputed boundary values in ψb to enforce proper boundary conditions in the discrete Laplacian operator.

Arguments

• ψb: Boundary buffer arrays for each component and direction.
• i: Index of the component being corrected.
• dir: Direction of the boundary (e.g., 1 = x, 2 = y, 3 = z).
• I: Cartesian index in the boundary array where the correction is computed.

Returns

A scalar value representing the diagonal Laplacian boundary correction at the given index.

multidomain_poisson!(ω, ψ, u, ψb, grid, fft_plan)

Solve Poisson equations across multiple grid levels using a multigrid-like approach.

Description

This function computes the solution of Poisson problems on a hierarchical set of grids. It performs the following steps for each level:

1. Coarsen the source term ω from finer to coarser grids.
2. Solve from coarse to fine.
3. Optionally compute the velocity field u at selected levels using curl!.

Arguments

• ω: Vector of vector fields representing source terms at each grid level.
• ψ: Vector of potential fields to store the solution at each level.
• u: Vector of velocity fields to update via curl of ψ.
• ψb: Boundary buffers for the potential fields.
• grid: Grid object containing the domain size and level information.
• fft_plan: Precomputed EigenbasisTransform plans for spectral solves.

Returns

The function updates ψ and u in-place with the Poisson solution and velocity field.

AbstractDeltaFunc

Abstract type for delta-function-like objects. Subtypes define specific delta kernels.

Usage

A delta function can be called on a vector r to evaluate the multidimensional delta:

delta(r)  # evaluates as the product of 1D delta values along each component

DeltaYang3S <: AbstractDeltaFunc

Smooth delta function approximation with compact support [-2, 2], following Yang et al. (2009).

• support(::DeltaYang3S) = 2 gives its support radius.
• Calling delta(r::Real) evaluates the function at a real point r using a piecewise formula:
  • |r| < 1 → first formula
  • 1 ≤ |r| < 2 → second formula
  • |r| ≥ 2 → returns 0

The function is smooth and satisfies partition-of-unity and moment conditions.

DeltaYang3S2 <: AbstractDeltaFunc

Smoother and wider delta function than DeltaYang3S, with compact support [-3, 3].

• support(::DeltaYang3S2) = 3 gives its support radius.
• Calling delta(x::Real) evaluates the function at x using piecewise formulas:
  • r ≤ 1
  • 1 < r ≤ 2
  • 2 < r ≤ 3
  • r > 3 → 0

Each segment uses polynomials, square roots, and arcsine terms to ensure smoothness and correct moment conditions.

Reg{D,T,N,A,M,W}
Reg(backend, T, delta, nb, Val{N})

Represents a regularization operator used for interpolation and spreading based on a discrete delta function.

Fields

• delta — the regularized delta function (a subtype of AbstractDeltaFunc).
• I — a matrix of index offsets defining the discrete stencil.
• weights — preallocated delta weights for each stencil point.

The struct is adapted for GPU execution via Adapt.@adapt_structure, allowing Reg objects to be transferred automatically between CPU and GPU memory.

Constructor

Reg(backend, T, delta, nb, Val{N}) creates a regularization operator in N dimensions. It allocates:

• the stencil index matrix I, and
• the multidimensional weights array whose size is determined by the support of the delta function.

backend controls where arrays are allocated (CPU or GPU), and nb is the number of bodies or markers for which weights are stored.

This type is typically used in immersed-boundary methods for evaluating and applying discrete delta functions.

update_weights!(reg, grid, xbs, ibs)

Update interpolation/spreading weights for immersed boundary markers.

This function computes the stencil indices and delta-function weights used to transfer data between Lagrangian marker positions (xbs) and the Eulerian grid (grid). Only markers listed in ibs are updated. The result is stored in-place inside the Reg object reg.

Arguments

• reg::Reg: Regularization structure containing delta kernel, stencil offsets, and a weight array to be filled.
• grid::Grid{N}: Eulerian grid used for mapping marker positions to grid coordinates.
• xbs: Array of marker positions (typically SVector{N,Float}).
• ibs: Indices of the markers to update.

Notes

• If ibs is empty, the function returns reg unchanged.
• For each marker and each velocity/force component, the function:
  1. Computes the integer grid offset I nearest to the marker.
  2. Iterates over all stencil points within the delta kernel’s support.
  3. Evaluates the delta function at normalized offsets (xb - xu) / h.
  4. Stores the resulting weights in reg.weights.

Returns

Returns the updated reg.

interpolate_body!(ub, reg, u)

Interpolate the Eulerian velocity field u onto the Lagrangian marker velocities ub using precomputed regularization stencils stored in reg.

This function gathers velocity values from the Eulerian grid for each marker and each velocity component, applies the corresponding delta–function weights, and stores the resulting interpolated velocities in-place in ub.

Arguments

• ub: Output array of marker velocities (e.g., Vector{SVector{N,T}}).
• reg::Reg: Regularization structure containing interpolation indices I and delta weights weights.
• u: Eulerian velocity field, given as an array of N grid arrays (u[1], u[2], …).

Notes

• For each marker, the function loops over velocity components and computes a weighted sum of nearby grid values using the delta kernel’s support.
• Uses the precomputed stencil offsets reg.I and weight tensors reg.weights, which must be updated before calling this function.
• Updates ub in-place and also returns it.

Returns

Updates ub in-place.

regularize!(fu, reg, fb)

Spread Lagrangian forces fb onto the Eulerian force field fu using the regularization stencils stored in reg.

This function distributes each marker force to nearby Eulerian grid points using the delta–function weights in reg.weights and the corresponding index offsets in reg.I. The resulting Eulerian force field is written in-place in fu.

Arguments

• fu: Output Eulerian force field, given as an array of N grids (fu[1], fu[2], …). All entries are reset to zero before accumulation.
• reg::Reg: Regularization structure containing interpolation/spreading indices I and delta weights weights.
• fb: Lagrangian forces, typically stored as Vector{SVector{N}}, one force vector per marker.

Notes

• For each marker, the force components are distributed over the delta kernel’s support region.
• This operation is the adjoint (transpose) of interpolate_body! in the immersed boundary method.
• Updates fu in-place and also returns it.

Returns

Updates fu in-place.