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:

Parameter	Type	Default	Range	Description
`profile`	enum	`"digamma"`	see above	Analytic profile function
`delta`	double	`1e-6`	>0 (floored)	Symmetric offset scale in response factor $(2\delta)$
`lambda`	double	`0.0`	unbounded	Scaling applied to warp response
`bias`	double	`0.0`	unbounded	Additive bias before profile evaluation

Power Profile:

Parameter	Type	Default	Range	Description
`exponent`	double	`2.0`	≥0	Power exponent for `profile = "power"`

Hyperexp/Qhyperexp:

Parameter	Type	Default	Range	Description
`hyperexp_epsilon`	double	`context epsilon` (fallback `1.0`)	>0	Pole shift control parameter
`hyperexp_depth`	integer	`context truncation level` (fallback `8`)	[1, 8192]	Truncation depth for sum
`hyperexp_q`	double	`0.9`	(0, 1)	q-deformation parameter

Complex Field Handling:

Parameter	Type	Default	Description
`complex_mode`	enum	`"component"`	`component`: process re/im independently; `polar`: preserve phase direction

Continuity Guards:

Parameter	Type	Default	Range	Description
`continuity_mode`	enum	`"none"`	`none`, `strict`, `clamped`, `limited`	Singularity handling policy
`continuity_clamp_min`	double	`-1e6`	unbounded	Lower bound for clamped mode
`continuity_clamp_max`	double	`1e6`	unbounded	Upper bound for clamped mode
`continuity_tolerance`	double	`1e-6`	>0 (floored)	Probe radius for guarded continuity modes

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

Parameter	Type	Default	Range	Description
`method`	enum	`"central"`	`central`, `forward`, `backward`	Finite difference stencil
`axis`	integer	`0`	[0, 1]	Derivative axis (0 = x, 1 = y)
`skew_forward`	boolean	`false`	—	Bias central scheme toward upwind
`spacing`	double	`1.0`	[1e-9, 10.0]	Grid spacing $\Delta x$ (alias: `dx`)
`boundary`	enum	`"periodic"`	see below	Boundary handling policy
`accumulate`	boolean	`false`	—	Add to output instead of overwriting

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

Parameter	Type	Default	Range	Description
`spacing_x`	double	`1.0`	[1e-9, ∞)	Grid spacing in X direction
`spacing_y`	double	`spacing_x`	[1e-9, ∞)	Grid spacing in Y direction (defaults to `spacing_x` if unset)
`axis_x`	integer	`-1` (auto)	[-1, 7]	Axis index for the X derivative (`-1` selects automatically)
`axis_y`	integer	`-1` (auto)	[-1, 7]	Axis index for the Y derivative (`-1` selects automatically; must differ from `axis_x`)
`stencil`	enum	`"central2"`	see below	Finite difference stencil type
`boundary`	enum	`"periodic"`	see below	Boundary condition handling
`accumulate`	boolean	`false`	—	Add to output instead of overwriting
`scale_by_dt`	boolean	`true`	—	Scale accumulated outputs by dt (`accumulate = true` only)

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

Parameter	Type	Default	Range	Description
`spacing_x`	double	`1.0`	[1e-9, ∞)	Grid spacing for $\partial v_x / \partial x$
`spacing_y`	double	`spacing_x`	[1e-9, ∞)	Grid spacing for $\partial v_y / \partial y$
`axis_x`	integer	`-1` (auto)	[-1, 7]	Axis index used for the X-derivative
`axis_y`	integer	`-1` (auto)	[-1, 7]	Axis index used for the Y-derivative
`stencil`	enum	`"central2"`	see below	Finite-difference stencil for both partials
`boundary`	enum	`"periodic"`	see below	Boundary sampling policy
`accumulate`	boolean	`false`	—	Add into the output instead of overwriting
`scale_by_dt`	boolean	`true`	—	Scale accumulated writes by dt (`accumulate = true` only)

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

Parameter	Type	Default	Range	Description
`spacing_x`	double	`1.0`	[1e-9, ∞)	Grid spacing for $\partial v_y / \partial x$
`spacing_y`	double	`spacing_x`	[1e-9, ∞)	Grid spacing for $\partial v_x / \partial y$
`axis_x`	integer	`-1` (auto)	[-1, 7]	Axis used for the X-derivative
`axis_y`	integer	`-1` (auto)	[-1, 7]	Axis used for the Y-derivative
`stencil`	enum	`"central2"`	see below	Finite-difference stencil for both partials
`boundary`	enum	`"periodic"`	see below	Boundary sampling policy
`accumulate`	boolean	`false`	—	Add into the output instead of overwriting
`scale_by_dt`	boolean	`true`	—	Scale accumulated writes by dt (`accumulate = true` only)

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
})