Advection Operators

🌀 Analytic Warp

sim_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

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

Returns: Operator handle (userdata)

Mathematical Formulation

The operator applies a profile function $\phi$ with symmetric gradient estimation:

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

where $\phi'$ is approximated via symmetric finite differences using the delta parameter.

Available profiles:

• digamma: $\phi(z) = \psi(z) = \frac{d}{dz}\ln\Gamma(z)$ — logarithmic derivative of gamma function
• trigamma: $\phi(z) = \psi_1(z) = \frac{d^2}{dz^2}\ln\Gamma(z)$ — second derivative
• power: $\phi(z) = z^p$ where $p$ is the exponent parameter
• tanh: $\phi(z) = \tanh(z)$ — hyperbolic tangent saturation
• hyperexp: $\phi(z) = \sum_{k=0}^{K-1} \frac{\lambda}{\lambda + \varepsilon + k}$ — hyperexponential sum
• qhyperexp: $\phi_q(z) = \sum_{k=0}^{K-1} \frac{\lambda}{\lambda + \varepsilon q^k}$ — q-deformed hyperexponential

Parameters

Core Parameters:

Parameter	Type	Default	Range	Description
profile	enum	"digamma"	see above	Analytic profile function
delta	double	1e-3	≥0	Symmetric offset for gradient estimation
lambda	double	1.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	1.0	≥0	Pole shift control parameter
hyperexp_depth	integer	8	≥1	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-9	≥0	Blend radius for limited mode

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
• The continuity_mode options help guard against singularities:
  • strict: Error on singular inputs
  • clamped: Clamp outputs to [continuity_clamp_min, continuity_clamp_max]
  • limited: Smooth blend near singularities using tolerance

Performance Notes

• Pointwise operation; scales linearly with field size
• hyperexp and qhyperexp profiles require $O(K)$ operations per sample where $K$ = hyperexp_depth
• Complex fields in polar mode require additional magnitude/phase conversions

Examples

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

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

-- Digamma warp with continuity protection
ooc.sim_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.sim_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.sim_add_analytic_warp_operator(ctx, field, {
    profile = "trigamma",
    delta = 1e-4,
    lambda = 0.05
})

---

∇ Spatial Derivative

sim_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

sim_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, the central scheme is biased toward the forward direction (useful for upwind schemes).

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 sim_add_gradient_operator instead

Runtime Configuration

Read current configuration:

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

Update configuration:

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

Examples

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

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

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

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

-- Query and modify at runtime
local cfg = ooc.sim_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.sim_spatial_derivative_update(ctx, op_index, {
    method = "backward",
    skew_forward = false,
    boundary = "neumann"
})

---

∇² Gradient

sim_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

sim_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}$
• 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	0	[0, 7]	Axis index for X derivative
axis_y	integer	1	[0, 7]	Axis index for Y derivative (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 writes by dt

Stencil options: forward2, backward2, central2, forward4, 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

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.sim_add_field(ctx, {256, 256}, { fill = 0.0 })
local uy = ooc.sim_add_field(ctx, {256, 256}, { fill = 0.0 })
ooc.sim_add_gradient_operator(ctx, u, ux, uy, {
    spacing_x = 0.1,
    spacing_y = 0.1
})

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

-- Upwind gradient for advection (forward in X, backward in Y)
ooc.sim_add_gradient_operator(ctx, u, ux, uy, {
    stencil = "forward2",
    boundary = "periodic"
})