Advection Operators

Analytic Warp 🌀

add_analytic_warp_operator(ctx, field, opts)

Apply nonlinear analytic deformations to field values using special function profiles. Useful for implementing nonlinear response curves, soft saturation, and phase-space warping.

Method Signature

add_analytic_warp_operator(ctx, field, [options]) -> operator

Returns: Operator handle (userdata)

Mathematical Formulation

For each sample, the operator computes a profile-gradient response and applies an explicit Euler update:

$$u \leftarrow u + \Delta t \cdot \lambda \cdot (2\delta) \cdot g(u + \text{bias})$$

where $g(\cdot)$ depends on profile:

• digamma: $g(z) = \psi_1(z)$ (trigamma)
• trigamma: $g(z) = \psi_2(z)$ (tetragamma)
• power: $g(z) = p|z|^{p-1}$
• tanh: $g(z) = 1 - \tanh^2(z)$
• hyperexp: $g(z) = \frac{d}{dz}\phi(z)$ with $\phi(z)=\sum_{k=0}^{K-1}\frac{z}{z+\varepsilon+k}$
• qhyperexp: $g(z)=\frac{d}{dz}\phi_q(z)$ with $\phi_q(z)=\sum_{k=0}^{K-1}\frac{z}{z+\varepsilon q^k}$

In continuity_mode = "clamped" or "limited", near singularities the operator may probe nearby samples to build a guarded response.

Parameters

Core Parameters:

Power Profile:

Hyperexp/Qhyperexp:

Complex Field Handling:

Continuity Guards:

Boundary & Initial Conditions

• Operates pointwise; no spatial boundary handling required
• Field values should be in a domain where the profile function is well-defined
• For digamma/trigamma, avoid non-positive integers (poles of gamma function)

Stability & Convergence

• Stability depends on the profile and lambda parameter
• Large lambda values can cause rapid field evolution; reduce timestep accordingly
• Lua-side default lambda = 0.0 means no evolution unless you set lambda explicitly
• The continuity_mode options help guard against singularities:
  • strict: Preserve direct analytic evaluation (no continuity fallback blend)
  • clamped: Clamp gradient responses to [continuity_clamp_min, continuity_clamp_max]
  • limited: Use guarded nearby probes and clamp within bounds

Performance Notes

• Pointwise operation; scales linearly with field size
• qhyperexp evaluation scales with truncation depth hyperexp_depth
• hyperexp uses digamma/trigamma differences and is typically less sensitive to hyperexp_depth
• Complex fields in polar mode require additional magnitude/phase conversions

Examples

-- Power-law nonlinearity with exponent 1.5
ooc.add_analytic_warp_operator(ctx, field, {
    profile = "power",
    exponent = 1.5,
    lambda = 0.4
})

-- Soft saturation using tanh
ooc.add_analytic_warp_operator(ctx, field, {
    profile = "tanh",
    lambda = 2.0,
    bias = 0.0
})

-- Digamma warp with continuity protection
ooc.add_analytic_warp_operator(ctx, field, {
    profile = "digamma",
    lambda = 0.1,
    continuity_mode = "clamped",
    continuity_clamp_min = -10.0,
    continuity_clamp_max = 10.0
})

-- Hyperexponential with q-deformation (complex field, polar mode)
ooc.add_analytic_warp_operator(ctx, complex_field, {
    profile = "qhyperexp",
    hyperexp_depth = 16,
    hyperexp_epsilon = 0.25,
    hyperexp_q = 0.85,
    complex_mode = "polar"
})

-- Trigamma profile for second-derivative response
ooc.add_analytic_warp_operator(ctx, field, {
    profile = "trigamma",
    delta = 1e-4,
    lambda = 0.05
})

---

Spatial Derivative ∇

add_spatial_derivative_operator(ctx, src, dst, opts)

Compute first-order spatial derivatives $\partial u / \partial x$ using finite difference stencils. Fundamental building block for advection, gradient computation, and PDE discretization.

Method Signature

add_spatial_derivative_operator(ctx, input, output, [options]) -> operator

Returns: Operator handle (userdata)

Mathematical Formulation

The operator computes $\frac{\partial u}{\partial x}$ (or $\frac{\partial u}{\partial y}$ when axis = 1) using finite differences:

Central difference (default):

$$\frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_{i+1} - u_{i-1}}{2\Delta x}$$

Forward difference:

$$\frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_{i+1} - u_i}{\Delta x}$$

Backward difference:

$$\frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_i - u_{i-1}}{\Delta x}$$

When skew_forward = true and method = "central", the runtime uses the forward stencil.

Parameters

Boundary options: periodic, neumann, dirichlet, reflective

Note: Lua accepts derivative as an alias for method.

Boundary & Initial Conditions

• Periodic: Wraps indices modulo field length
• Neumann: Zero normal derivative; ghost cells mirror interior values
• Dirichlet: Zero values at boundaries
• Reflective: Mirror-reflects values across boundary

Stability & Convergence

• Pure derivative operators are unconditionally stable
• Truncation error:
  • Central: $O(\Delta x^2)$
  • Forward/Backward: $O(\Delta x)$
• For advection problems, use upwind schemes (forward for positive velocity, backward for negative) to ensure stability
• Central schemes can introduce dispersion errors in advection; skew_forward provides a compromise

Performance Notes

• Nonlocal operation requiring neighbor access
• Does not use dt scaling by default (pure spatial operator)
• For 2D gradients requiring both $\partial/\partial x$ and $\partial/\partial y$, consider add_gradient_operator instead

Runtime Configuration

Read current configuration:

local cfg = ooc.spatial_derivative_config(ctx, op_index)
-- Returns: { input_field, output_field, spacing, method, method_index, axis, skew_forward, boundary, accumulate }

Update configuration:

ooc.spatial_derivative_update(ctx, op_index, {
    axis = 1,
    method = "backward",
    spacing = 0.05,
    boundary = "neumann"
})

Examples

-- Basic central difference
ooc.add_spatial_derivative_operator(ctx, u, du_dx, {
    method = "central",
    spacing = 0.1
})

-- Forward difference for upwind advection (positive velocity)
ooc.add_spatial_derivative_operator(ctx, u, du_dx, {
    method = "forward",
    spacing = 0.05,
    boundary = "periodic"
})

-- Y-derivative with Dirichlet boundaries
ooc.add_spatial_derivative_operator(ctx, u, du_dy, {
    method = "central",
    axis = 1,
    spacing = 0.1,
    boundary = "dirichlet"
})

-- Accumulate derivative into existing field
ooc.add_spatial_derivative_operator(ctx, u, rhs, {
    method = "backward",
    spacing = dx,
    accumulate = true
})

-- Query and modify at runtime
local cfg = ooc.spatial_derivative_config(ctx, op_index)
if cfg then
    ooc.log("∂/∂x: axis=%d method=%s dx=%.3g", cfg.axis, cfg.method, cfg.spacing)
end

ooc.spatial_derivative_update(ctx, op_index, {
    method = "backward",
    skew_forward = false,
    boundary = "neumann"
})

---

Gradient ∇

add_gradient_operator(ctx, input, output_x, output_y, opts)

Compute 2D gradient components $(\partial_x u, \partial_y u)$ using finite difference stencils. Outputs the X and Y partial derivatives to separate fields.

Method Signature

add_gradient_operator(ctx, input, output_x, output_y, [options]) -> operator

Returns: Operator handle (userdata)

Mathematical Formulation

For a scalar field $u(x, y)$, the gradient is:

$$\nabla u = \left( \frac{\partial u}{\partial x}, \frac{\partial u}{\partial y} \right)$$

The discrete approximation depends on the selected stencil:

• central2 (2nd-order central): $\frac{\partial u}{\partial x} \approx \frac{u_{i+1} - u_{i-1}}{2 \Delta x}$
• central4 (4th-order central): $\frac{\partial u}{\partial x} \approx \frac{-u_{i+2} + 8u_{i+1} - 8u_{i-1} + u_{i-2}}{12 \Delta x}$
• forward1 (1st-order forward): $\frac{\partial u}{\partial x} \approx \frac{u_{i+1} - u_i}{\Delta x}$
• backward1 (1st-order backward): $\frac{\partial u}{\partial x} \approx \frac{u_i - u_{i-1}}{\Delta x}$
• forward2 (2nd-order forward): $\frac{\partial u}{\partial x} \approx \frac{-3u_i + 4u_{i+1} - u_{i+2}}{2 \Delta x}$
• backward2 (2nd-order backward): $\frac{\partial u}{\partial x} \approx \frac{3u_i - 4u_{i-1} + u_{i-2}}{2 \Delta x}$

Parameters

Stencil options: forward1, backward1, forward2, backward2, central2, central4

Boundary options: periodic, neumann, dirichlet, reflective

Boundary & Initial Conditions

• Periodic: Wraps around domain edges using modular indexing
• Neumann: Zero-gradient at boundaries (ghost cells mirror interior)
• Dirichlet: Zero values at boundaries
• Reflective: Mirror-reflects values at boundaries

Stability & Convergence

• Central difference schemes are unconditionally stable for pure gradient computation
• Higher-order stencils (central4) reduce truncation error from $O(\Delta x^2)$ to $O(\Delta x^4)$
• Forward/backward stencils introduce directional bias; use for upwind schemes in advection problems
• A single stencil setting is used for both gradient components

Performance Notes

• Outputs to two separate fields; ensure both are allocated before calling
• Complex fields are processed component-wise (real and imaginary gradients computed independently)
• Stencil width affects memory access patterns; central4 requires wider halos

Examples

-- Basic 2D gradient with default central differences
local ux = ooc.add_field(ctx, {256, 256}, { fill = 0.0 })
local uy = ooc.add_field(ctx, {256, 256}, { fill = 0.0 })
ooc.add_gradient_operator(ctx, u, ux, uy, {
    spacing_x = 0.1,
    spacing_y = 0.1
})

-- High-order gradient with Neumann boundaries
ooc.add_gradient_operator(ctx, u, ux, uy, {
    stencil = "central4",
    spacing_x = 0.05,
    boundary = "neumann"
})

-- Upwind-style gradient for advection (single stencil applies to both axes)
ooc.add_gradient_operator(ctx, u, ux, uy, {
    stencil = "forward2",
    boundary = "periodic"
})

---

Divergence ∇·

add_divergence_operator(ctx, vx, vy, out, opts)

Compute the scalar divergence of a 2D vector field using matched finite-difference stencils on the X and Y components.

Method Signature

add_divergence_operator(ctx, input_x, input_y, output, [options]) -> operator

Returns: Operator handle (userdata)

Mathematical Formulation

For a vector field $\mathbf{v} = (v_x, v_y)$:

$$\nabla \cdot \mathbf{v} = \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y}$$

The selected stencil is applied independently to each partial derivative before summation:

• central2: $\frac{f_{i+1} - f_{i-1}}{2\Delta}$
• central4: $\frac{-f_{i+2} + 8f_{i+1} - 8f_{i-1} + f_{i-2}}{12\Delta}$
• forward1/backward1: first-order one-sided derivatives
• forward2/backward2: second-order one-sided derivatives

Parameters

Stencil options: forward1, backward1, forward2, backward2, central2, central4

Boundary options: periodic, neumann, dirichlet, reflective

Boundary & Initial Conditions

• Uses the selected boundary policy when sampling neighbors
• axis_x and axis_y must resolve to different axes
• Complex fields are differentiated component-wise before summing

Stability & Convergence

• Pure divergence evaluation is unconditionally stable
• central4 improves truncation error over central2 for smooth fields
• One-sided stencils introduce directional bias and are best reserved for upwind-style workflows

Performance Notes

• Reads two input fields and writes one output field
• Wider stencils require broader halos and more neighbor fetches
• accumulate = true is useful when divergence is one term within a larger RHS assembly

Examples

-- Incompressibility diagnostic
ooc.add_divergence_operator(ctx, vx, vy, div, {
    stencil = "central4",
    spacing_x = 0.5,
    spacing_y = 0.5
})

-- Nonuniform grid spacing with Neumann boundaries
ooc.add_divergence_operator(ctx, flux_x, flux_y, rhs, {
    spacing_x = 0.1,
    spacing_y = 0.25,
    boundary = "neumann"
})

-- Accumulate divergence contribution into an existing RHS field
ooc.add_divergence_operator(ctx, vx, vy, rhs, {
    accumulate = true,
    scale_by_dt = true
})

---

Curl ∇×

add_curl_operator(ctx, vx, vy, out, opts)

Compute the scalar 2D curl (vorticity) of a vector field using the same stencil family as the gradient and divergence operators.

Method Signature

add_curl_operator(ctx, input_x, input_y, output, [options]) -> operator

Returns: Operator handle (userdata)

Mathematical Formulation

For a planar vector field $\mathbf{v} = (v_x, v_y)$, the scalar curl is:

$$(\nabla \times \mathbf{v})_z = \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y}$$

Each partial derivative uses the selected stencil and its corresponding spacing:

• spacing_x scales $\partial v_y / \partial x$
• spacing_y scales $\partial v_x / \partial y$

Parameters

Stencil options: forward1, backward1, forward2, backward2, central2, central4

Boundary options: periodic, neumann, dirichlet, reflective

Boundary & Initial Conditions

• Uses the selected boundary mode when sampling off-edge neighbors
• axis_x and axis_y must resolve to different axes
• Complex fields yield component-wise curl outputs

Stability & Convergence

• Curl evaluation itself is unconditionally stable
• Centered stencils minimize bias for vorticity diagnostics
• One-sided stencils are useful when matching a broader directional discretization scheme

Performance Notes

• Reads two component fields and writes one scalar output field
• Shares the same stencil cost profile as add_divergence_operator
• Useful for vorticity visualization, circulation estimates, and rotational forcing terms

Examples

-- Basic vorticity field
ooc.add_curl_operator(ctx, vx, vy, omega, {
    stencil = "central4",
    spacing_x = 0.5,
    spacing_y = 0.5
})

-- Periodic curl on a square lattice
ooc.add_curl_operator(ctx, flow_x, flow_y, vort, {
    boundary = "periodic",
    spacing_x = 0.1
})

-- Add rotational term into an existing buffer
ooc.add_curl_operator(ctx, vx, vy, rhs, {
    accumulate = true,
    scale_by_dt = true
})