Diffusion Operators

Linear Dissipative 📉

add_linear_dissipative_operator(ctx, field, opts)

Apply spectral fractional Laplacian diffusion via FFT. Real and complex fields share the same real spectral multiplier, so real inputs remain real while still benefiting from the spectral update.

Method Signature

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

Returns: Operator handle (userdata)

Mathematical Formulation

The operator applies damping in spectral space:

\hat{\psi}(k) \leftarrow \hat{\psi}(k) \cdot \exp\left(-\nu |k|^\alpha \Delta t\right)

where:

$\hat{\psi}(k)$ is the Fourier coefficient at wavenumber $k$
$\nu$ is the viscosity parameter
$\alpha$ is the fractional Laplacian order (2 = classical diffusion)
$|k|$ is scaled by the spacing parameter

Special cases:

$\alpha = 2$ : Classical diffusion (heat equation)
$\alpha = 1$ : Half-Laplacian / fractional diffusion
$\alpha > 2$ : Hyperdiffusive-style high-frequency damping
$0 < \alpha < 2$ : Fractional diffusion with heavier tails than the classical heat flow

Parameters

Parameter	Type	Default	Range	Description
`viscosity`	double	`0.0`	[0, 10]	Diffusion coefficient $\nu$ (required)
`alpha`	double	`2.0`	[0, 10]	Fractional Laplacian power
`spacing`	double	`1.0`	[1e-9, 10]	Grid spacing for wavenumber scaling

Boundary & Initial Conditions

Periodic boundaries only: Spectral methods assume periodic domain
For non-periodic boundaries, use add_laplacian_operator instead
Field should be initialized before first application
Real and complex fields are both supported

Stability & Convergence

Unconditionally stable for any timestep (exponential damping)
viscosity = 0.0 results in a no-op pass-through
Larger alpha values produce sharper high-frequency damping
Spectral accuracy in space; exponential integrator in time

Performance Notes

Requires an FFT/IFFT pair when starting from physical space (O(N log N))
Real inputs use the same real spectral multiplier; explicit complex storage is not required
Most efficient for periodic domains with power-of-two sizes
Consider fusing with other spectral operators using add_linear_spectral_fusion_operator

Examples

-- Classical diffusion (alpha = 2)
ooc.add_linear_dissipative_operator(ctx, field, {
    viscosity = 0.5,
    spacing = 0.1
})

-- Fractional diffusion for anomalous transport
ooc.add_linear_dissipative_operator(ctx, field, {
    viscosity = 0.2,
    alpha = 1.5,  -- superdiffusion
    spacing = 0.05
})

-- Hyperdiffusion for numerical stability
ooc.add_linear_dissipative_operator(ctx, field, {
    viscosity = 1e-4,
    alpha = 4.0  -- targets high-k modes
})

-- No-op passthrough
ooc.add_linear_dissipative_operator(ctx, field, {
    viscosity = 0.0
})

Linear Spectral Fusion 🔮

add_linear_spectral_fusion_operator(ctx, field, opts)

Fuse dissipation, dispersion, and phase rotation in a single spectral step. More efficient than chaining separate spectral operators when multiple effects are needed.

Method Signature

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

Returns: Operator handle (userdata)

Mathematical Formulation

The operator combines three spectral effects:

\hat{\psi}(k) \leftarrow \hat{\psi}(k) \cdot \exp\left(-\nu|k|^\alpha \Delta t\right) \cdot \exp\left(i \left(\beta \left||k| - k_0\right|^p \Delta t + \Omega_\phi \Delta t\right)\right)

Components:

Dissipation: $\exp(-\nu|k|^\alpha \Delta t)$ — amplitude damping
Dispersion: $\exp(i \beta ||k| - k_0|^p \Delta t)$ — wavenumber-dependent phase velocity
Phase rotation: $\exp(i \Omega_\phi \Delta t)$ — uniform complex rotation

For real-valued fields, the runtime uses the projected real multiplier

\exp\left(-\nu |k|^\alpha \Delta t\right)\cos\left(\beta \left||k| - k_0\right|^p \Delta t + \Omega_\phi \Delta t\right)

so the result stays in the real subspace.

Parameters

Dissipation:

Parameter	Type	Default	Range	Description
`viscosity`	double	`0.0`	≥0	Dissipation strength $\nu$
`alpha`	double	`2.0`	[0, 2]	Fractional Laplacian exponent
`dissipation_spacing`	double	`1.0`	[1e-9, 1]	Grid spacing for dissipation term

Dispersion:

Parameter	Type	Default	Range	Description
`dispersion_coefficient`	double	`0.0`	[-1, 1]	Phase-rate amplitude $\beta$
`dispersion_order`	double	`1.0`	[0, 4]	Power $p$ on $
`dispersion_reference_k`	double	`0.0`	[0, 1]	Reference wavenumber $k_0$
`dispersion_spacing`	double	`1.0`	[1e-9, 1]	Grid spacing for dispersion

Phase Rotation:

Parameter	Type	Default	Range	Description
`phase_rate`	double	`0.0`	[0, 10]	Global phase rate $\Omega_\phi$ (rad/s)

Boundary & Initial Conditions

Periodic boundaries required (spectral method)
Real and complex fields are both supported
All three effects can be enabled/disabled independently by setting coefficients to zero

Stability & Convergence

Complex fields combine unitary phase rotation with exponential damping
Real fields use the projected cosine factor and remain in the real subspace
Dispersion with dispersion_order = 2 models group velocity dispersion
dispersion_order = 3 adds third-order dispersion (pulse asymmetry)

Performance Notes

Single FFT/IFFT pair when starting from physical space, regardless of how many effects are enabled
More efficient than chaining linear_dissipative + dispersion + phase_rotate
Real fields use the projected real multiplier; complex fields use the full complex phase factor
Recommended for systems requiring multiple spectral operations

Examples

-- Combined dissipation and dispersion (Schrödinger-like)
ooc.add_linear_spectral_fusion_operator(ctx, field, {
    viscosity = 0.1,
    alpha = 2.0,
    dispersion_coefficient = 1.0,
    dispersion_order = 2.0
})

-- Pure dispersion with carrier offset
ooc.add_linear_spectral_fusion_operator(ctx, field, {
    dispersion_coefficient = 0.5,
    dispersion_order = 1.0,
    dispersion_reference_k = 0.3
})

-- Dissipation + global phase rotation (detuned oscillator)
ooc.add_linear_spectral_fusion_operator(ctx, field, {
    viscosity = 0.05,
    phase_rate = 2.0 * math.pi  -- 1 Hz rotation
})

-- Full fusion: dissipation + dispersion + rotation
ooc.add_linear_spectral_fusion_operator(ctx, field, {
    viscosity = 0.2,
    alpha = 2.0,
    dissipation_spacing = 0.1,
    dispersion_coefficient = 1.0,
    dispersion_order = 2.0,
    dispersion_reference_k = 0.5,
    phase_rate = 0.1
})

Dispersion 🌈

add_dispersion_operator(ctx, field, opts)

Apply pure spectral phase modulation for modeling wave dispersion. No amplitude change; only phase velocities are affected.

Method Signature

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

Returns: Operator handle (userdata)

Mathematical Formulation

\hat{\psi}(k) \leftarrow \hat{\psi}(k) \cdot \exp\left(i \beta \left||k| - k_0\right|^p \Delta t\right)

where:

$\beta$ is the coefficient (can be positive or negative)
$p$ is the order parameter
$k_0$ is the reference_k (carrier wavenumber offset)

Physical interpretation:

order = 1: Linear dispersion (constant group velocity offset)
order = 2: Group velocity dispersion (pulse broadening)
order = 3: Third-order dispersion (pulse asymmetry)
Negative coefficient reverses dispersion direction

Parameters

Parameter	Type	Default	Range	Description
`coefficient`	double	`0.0`	[-10, 10]	Dispersion strength $\beta$ (required)
`order`	double	`2.0`	[0, 10]	Power on $
`spacing`	double	`1.0`	[1e-9, 1]	Grid spacing for wavenumber
`reference_k`	double	`0.0`	[0, 1]	Reference wavenumber $k_0$

Boundary & Initial Conditions

Periodic boundaries required
Real and complex fields are both supported
Real fields use the projected cosine multiplier instead of the full complex phase factor

Stability & Convergence

Complex fields: Pure phase rotation is unitary and energy-conserving
Real fields: The projected cosine factor keeps outputs real but can modulate amplitude
Higher order values create stronger k-dependent phase velocities
reference_k > 0 shifts the dispersion relation away from k=0

Performance Notes

Requires an FFT/IFFT pair when starting from physical space
Real inputs stay real by using the projected cosine multiplier
Consider add_linear_spectral_fusion_operator if combining with dissipation

Examples

-- Quadratic dispersion (group velocity dispersion)
ooc.add_dispersion_operator(ctx, field, {
    coefficient = 1.0,
    order = 2.0
})

-- Negative dispersion (anomalous regime)
ooc.add_dispersion_operator(ctx, field, {
    coefficient = -0.5,
    order = 2.0
})

-- Linear dispersion with carrier offset
ooc.add_dispersion_operator(ctx, field, {
    coefficient = 1.0,
    order = 1.0,
    reference_k = 0.25
})

-- Third-order dispersion for pulse shaping
ooc.add_dispersion_operator(ctx, field, {
    coefficient = 0.1,
    order = 3.0,
    spacing = 0.05
})

Fractional Memory 🧠

add_fractional_memory_operator(ctx, field, opts)

Apply long-range temporal memory using fractional calculus. Models systems with power-law relaxation, viscoelastic materials, and anomalous diffusion.

Method Signature

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

Returns: Operator handle (userdata)

Mathematical Formulation

The operator implements a Grünwald-Letnikov fractional derivative:

u^{n+1} = u^n + g \cdot \Delta t^{1-q} \sum_{j=0}^{M-1} c_j \cdot u^{n-j}

where:

$q$ is the fractional order ( $0 < q \leq 1$ )
$g$ is the gain parameter
$M$ is memory_steps
$c_j$ are the Grünwald-Letnikov coefficients: $c_j = (-1)^j \binom{q}{j}$

Physical interpretation:

$q = 1$ : Standard first-order derivative (exponential decay)
$q < 1$ : Sub-diffusive memory (power-law relaxation)
Lower $q$ = longer memory, slower decay

Parameters

Parameter	Type	Default	Range	Description
`order`	double	(0, 1]	Fractional derivative order $q$ (required)
`gain`	double	`0.5`	unbounded	Scale applied to fractional response
`memory_steps`	integer	`32`	≥1	Number of past samples retained

Boundary & Initial Conditions

Operates pointwise; no spatial boundaries
Maintains internal history buffer of memory_steps samples
First memory_steps timesteps may show transient behavior as history fills

Stability & Convergence

Stability depends on gain, order, and timestep
Smaller order values require more memory_steps for accuracy
Rule of thumb: memory_steps >= 10 / order for reasonable accuracy
Large gain values can cause instability; start small and increase

Performance Notes

Memory usage: O(memory_steps × field_size)
Computation: O(memory_steps) per sample per timestep
History is stored in a circular buffer; no reallocation during runtime
Consider reducing memory_steps for real-time applications

Runtime Configuration

Read current configuration:

local cfg = ooc.fractional_memory_config(ctx, op_index)
-- Returns: { field_index, order, gain, memory_steps }

Update configuration:

ooc.fractional_memory_update(ctx, op_index, {
    memory_steps = 256,
    gain = 0.35
})

Examples

-- Standard fractional memory (order 0.6)
ooc.add_fractional_memory_operator(ctx, field, {
    order = 0.6,
    gain = 0.25,
    memory_steps = 128
})

-- Near-classical first-order lag
ooc.add_fractional_memory_operator(ctx, field, {
    order = 1.0,
    gain = 0.5
})

-- Long-range memory for viscoelastic modeling
ooc.add_fractional_memory_operator(ctx, field, {
    order = 0.3,  -- strong memory effect
    gain = 0.1,
    memory_steps = 512
})

-- Runtime adjustment
local cfg = ooc.fractional_memory_config(ctx, op_index)
if cfg then
    ooc.log("fracmem: order=%.2f gain=%.3f steps=%d",
        cfg.order, cfg.gain, cfg.memory_steps)
end

ooc.fractional_memory_update(ctx, op_index, {
    memory_steps = 256,
    gain = 0.35
})

Laplacian ∇²

add_laplacian_operator(ctx, input, output, opts)

Compute the discrete Laplacian $\nabla^2 u$ using cross or isotropic stencils. Essential for diffusion equations, heat transfer, and Poisson problems.

Method Signature

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

Returns: Operator handle (userdata)

Mathematical Formulation

The Laplacian operator in 2D is:

\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}

The discrete approximation depends on the selected stencil:

Cross stencils (axis-aligned):

cross2 (5-point stencil):

\nabla^2 u \approx \frac{u_{i+1,j} + u_{i-1,j} - 2u_{i,j}}{\Delta x^2} + \frac{u_{i,j+1} + u_{i,j-1} - 2u_{i,j}}{\Delta y^2}

cross4 (9-point extended):

\nabla^2 u \approx \frac{-u_{i+2} + 16u_{i+1} - 30u_i + 16u_{i-1} - u_{i-2}}{12\Delta x^2} + \text{(y terms)}

Isotropic stencil (rotationally symmetric for $\Delta x = \Delta y$ ):

isotropic9 (9-point): Includes diagonal neighbors for better rotational symmetry:

\nabla^2 u \approx \frac{1}{6h^2}\left(u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} + u_{i+1,j+1} + u_{i-1,j-1} + u_{i+1,j-1} + u_{i-1,j+1} - 8u_{i,j}\right)

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 the X term (`-1` selects automatically)
`axis_y`	integer	`-1` (auto)	[-1, 7]	Axis index for the Y term (`-1` selects automatically; must differ from `axis_x`)
`stencil`	enum	`"cross2"`	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 (`accumulate = true` only)

Stencil options: cross2, cross4, isotropic9

Boundary options: periodic, neumann, dirichlet, reflective

Boundary & Initial Conditions

Periodic: Wraps around domain edges; suitable for toroidal topologies
Neumann: Zero normal derivative at boundaries (∂u/∂n = 0); models insulated boundaries
Dirichlet: Fixed values at boundaries (u = 0); models grounded or fixed-temperature edges
Reflective: Mirror boundary values; creates symmetric extension
axis_x and axis_y must resolve to different axes when set explicitly

Stability & Convergence

When used with explicit time integration (e.g., Euler), the CFL condition for diffusion is:

\Delta t \leq \frac{(\Delta x)^2}{4D}

where $D$ is the diffusion coefficient. Implicit methods (backward Euler, Crank-Nicolson) remove this constraint.

Truncation error:

cross2: $O(\Delta x^2)$
cross4: $O(\Delta x^4)$
isotropic9: $O(\Delta x^2)$ but with improved rotational symmetry

Performance Notes

isotropic9 requires equal spacing ( $\Delta x = \Delta y$ ); raises an error otherwise
cross4 requires larger halo regions, increasing memory bandwidth
For complex fields, real and imaginary parts are processed independently
Consider add_linear_dissipative_operator for spectral (FFT-based) diffusion when periodic boundaries suffice

Examples

-- Basic 2D Laplacian with 5-point stencil
local lap = ooc.add_field(ctx, {256, 256}, { fill = 0.0 })
ooc.add_laplacian_operator(ctx, u, lap, {
    spacing_x = 0.1,
    spacing_y = 0.1
})

-- High-order Laplacian for smooth fields
ooc.add_laplacian_operator(ctx, u, lap, {
    stencil = "cross4",
    spacing_x = 0.05,
    boundary = "neumann"
})

-- Isotropic 9-point stencil (requires equal spacing)
ooc.add_laplacian_operator(ctx, u, lap, {
    stencil = "isotropic9",
    spacing_x = 0.1,
    spacing_y = 0.1,  -- must equal spacing_x
    boundary = "dirichlet"
})

-- Heat equation: du/dt = D * laplacian(u)
local D = 0.01
ooc.add_laplacian_operator(ctx, temp, dtemp, {
    stencil = "cross2",
    spacing_x = dx,
    accumulate = true
})
ooc.add_scale_operator(ctx, dtemp, dtemp, { scale = D })