Diffusion Operators

📉 Linear Dissipative

sim_add_linear_dissipative_operator(ctx, field, opts)

Apply spectral fractional Laplacian diffusion via FFT. Provides efficient, unconditionally stable damping for periodic domains with configurable fractional order.

Method Signature

sim_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$: Hyperdiffusion
• $\alpha < 2$: Superdiffusion (Lévy flights)

Parameters

Parameter	Type	Default	Range	Description
viscosity	double	[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 sim_add_laplacian_operator instead
• Field should be initialized before first application

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 FFT/IFFT pair per application (O(N log N))
• Complex storage required; real fields are promoted internally
• Most efficient for periodic domains with power-of-two sizes
• Consider fusing with other spectral operators using sim_add_linear_spectral_fusion_operator

Examples

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

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

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

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

---

🔮 Linear Spectral Fusion

sim_add_linear_spectral_fusion_operator(ctx, field, opts)

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

Method Signature

sim_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

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, 10]	Grid spacing for dissipation term

Dispersion:

Parameter	Type	Default	Range	Description
dispersion_coefficient	double	0.0	[-10, 10]	Phase-rate amplitude $\beta$
dispersion_order	double	1.0	[0, 10]	Power $p$ on $
dispersion_reference_k	double	0.0	[0, 10]	Reference wavenumber $k_0$
dispersion_spacing	double	1.0	[1e-9, 10]	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)
• Complex field storage required
• All three effects can be enabled/disabled independently by setting coefficients to zero

Stability & Convergence

• Unconditionally stable (unitary phase rotation, exponential damping)
• Dispersion with dispersion_order = 2 models group velocity dispersion
• dispersion_order = 3 adds third-order dispersion (pulse asymmetry)

Performance Notes

• Single FFT/IFFT pair regardless of how many effects are enabled
• More efficient than chaining linear_dissipative + dispersion + phase_rotate
• Recommended for systems requiring multiple spectral operations

Examples

-- Combined dissipation and dispersion (Schrödinger-like)
ooc.sim_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.sim_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.sim_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.sim_add_linear_spectral_fusion_operator(ctx, field, {
    viscosity = 0.2,
    alpha = 2.0,
    dissipation_spacing = 0.1,
    dispersion_coefficient = 1.5,
    dispersion_order = 2.0,
    dispersion_reference_k = 0.5,
    phase_rate = 0.1
})

---

🌈 Dispersion

sim_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

sim_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	[-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
• Complex field storage required
• Preserves total field energy (unitary operation)

Stability & Convergence

• Unconditionally stable (pure phase rotation is unitary)
• Energy-conserving; no amplitude damping
• Higher order values create stronger k-dependent phase velocities
• reference_k > 0 shifts the dispersion relation away from k=0

Performance Notes

• Requires FFT/IFFT pair
• Consider sim_add_linear_spectral_fusion_operator if combining with dissipation

Examples

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

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

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

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

---

🧠 Fractional Memory

sim_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

sim_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.sim_fractional_memory_config(ctx, op_index)
-- Returns: { field_index, order, gain, memory_steps }

Update configuration:

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

Examples

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

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

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

-- Runtime adjustment
local cfg = ooc.sim_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.sim_fractional_memory_update(ctx, op_index, {
    memory_steps = 256,
    gain = 0.35
})

---

∇² Laplacian

sim_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

sim_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 stencils (rotationally symmetric for $\Delta x = \Delta y$):

• isotropic5 (5-point): Standard cross stencil with equal weights
• 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	0	[0, 7]	Axis index for X term
axis_y	integer	1	[0, 7]	Axis index for Y term (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

Stencil options: cross2, cross4, isotropic5, 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

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 sim_add_linear_dissipative_operator for spectral (FFT-based) diffusion when periodic boundaries suffice

Examples

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

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

-- Isotropic 9-point stencil (requires equal spacing)
ooc.sim_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.sim_add_laplacian_operator(ctx, temp, dtemp, {
    stencil = "cross2",
    spacing_x = dx,
    accumulate = true
})
ooc.sim_add_scale_operator(ctx, dtemp, dtemp, { scale = D })