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, 10.0]	Grid spacing in X direction
spacing_y	double	spacing_x	[1e-9, 10.0]	Grid spacing in Y direction (defaults to spacing_x if unset)
axis_x	integer	-1 (auto)	[-1, 7]	Axis index for X derivative (-1 selects last axis)
axis_y	integer	-1 (auto)	[-1, 7]	Axis index for Y derivative (-1 selects second-last axis; 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 outputs by dt

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