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

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:

Dispersion:

Phase Rotation:

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

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

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 $\Delta u = \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,j} + 16u_{i+1,j} - 30u_{i,j} + 16u_{i-1,j} - u_{i-2,j}}{12\Delta x^2} + \frac{-u_{i,j+2} + 16u_{i,j+1} - 30u_{i,j} + 16u_{i,j-1} - u_{i,j-2}}{12\Delta y^2}$$

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(4\left(u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1}\right) + u_{i+1,j+1} + u_{i-1,j-1} + u_{i+1,j-1} + u_{i-1,j+1} - 20u_{i,j}\right)$$

Parameters

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

For the standard 2D cross2 stencil under explicit time integration (e.g., Euler), the diffusion CFL condition is:

$$\Delta t \leq \frac{1}{2D\left(\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}\right)}$$

For equal spacing $\Delta x = \Delta y = h$, this reduces to

$$\Delta t \leq \frac{h^2}{4D}$$

where $D$ is the diffusion coefficient. Higher-order or isotropic stencils can have different stability limits. 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 })