Custom Operators

Creating custom differential and integral operators.

Overview

Beyond built-in operators (grad, div, curl, lap), you can define custom operators for specialized physics.

Helper Functions

Simple Custom Operators

# Advection operator: u·∇f
function advection(u, f)
    result = u.components[1].data .* ∂x(f).data
    for i in 2:length(u.components)
        result .+= u.components[i].data .* d_operators[i](f).data
    end
    return result
end

Using in Equations

# Option 1: Expand manually
Tarang.add_equation!(problem, "∂t(T) = -ux*∂x(T) - uz*∂z(T)")

# Option 2: Use helper in equation string
# (requires registration)

Using Operators in Equations

The equation parser recognizes all built-in operators. Use them directly:

# Built-in operators available in equations:
# grad, div, curl, lap (or Δ), dt (or ∂t), d
# integrate, average, interpolate, convert, lift
# sin, cos, tan, exp, log, sqrt, abs, tanh

# Use operators directly in equations
Tarang.add_equation!(problem, "∂t(T) - kappa*Δ(T) = -ux*∂x(T) - uz*∂z(T)")

Spectral Differentiation

Fourier Derivative

function fourier_derivative(field, basis, order=1)
    k = get_wavenumbers(basis)
    Tarang.ensure_layout!(field, :c)

    # Multiply by (ik)^order
    get_coeff_data(field) .*= (1im .* k) .^ order

    return field
end

Chebyshev Derivative

function chebyshev_derivative(field, basis)
    Tarang.ensure_layout!(field, :c)

    # Use differentiation matrix
    D = chebyshev_diff_matrix(basis.size)
    get_coeff_data(field) .= D * get_coeff_data(field)

    return field
end

Integral Operators

Spatial Integration

function integrate(field, dim)
    Tarang.ensure_layout!(field, :g)

    if dim == 1  # x-direction
        dx = field.bases[1].length / field.bases[1].size
        return sum(get_grid_data(field), dims=1) * dx
    elseif dim == 2  # z-direction
        # Chebyshev: use quadrature weights
        weights = chebyshev_weights(field.bases[2])
        return sum(get_grid_data(field) .* weights', dims=2)
    end
end

Running Average

function running_average(field, window)
    Tarang.ensure_layout!(field, :g)

    # Convolution with box filter
    # ...
end

Vector Calculus

Strain Rate Tensor

function strain_rate(u)
    # S_ij = 1/2 (∂u_i/∂x_j + ∂u_j/∂x_i)
    S = TensorField(u.dist, u.coords, "S", u.bases; symmetric=true)

    ux, uz = u.components[1], u.components[2]

    S[1,1] = ∂x(ux)
    S[1,2] = 0.5 * (∂x(uz) + ∂z(ux))
    S[2,2] = ∂z(uz)

    return S
end

Vorticity (2D)

function vorticity_2d(u)
    ux, uz = u.components[1], u.components[2]
    return ∂x(uz) - ∂z(ux)
end

Helicity (3D)

function helicity(u)
    omega = ∇×(u)

    # H = u · ω
    H = u.components[1].data .* omega.components[1].data
    for i in 2:3
        H .+= u.components[i].data .* omega.components[i].data
    end

    return H
end

Nonlinear Operators

Convective Derivative

function convective_derivative(u, f)
    # (u·∇)f
    result = zeros(size(get_grid_data(f)))

    for (i, comp) in enumerate(u.components)
        Tarang.ensure_layout!(comp, :g)
        df = d_operators[i](f)
        Tarang.ensure_layout!(df, :g)
        result .+= get_grid_data(comp) .* get_grid_data(df)
    end

    return result
end

Nonlinear Term (u·∇u)

function nonlinear_advection(u)
    # Returns vector field
    result = similar(u)

    for (j, uj) in enumerate(u.components)
        get_grid_data(result.components[j]) .= convective_derivative(u, uj)
    end

    return result
end

Physics-Specific Operators

Coriolis Force

function coriolis(u, Omega)
    # 2Ω × u
    # For rotation about z-axis:
    fx = -2 * Omega * u.components[2]  # -2Ω*v
    fy =  2 * Omega * u.components[1]  # +2Ω*u
    fz = 0

    return (fx, fy, fz)
end

Lorentz Force (MHD)

function lorentz_force(J, B)
    # J × B
    return cross(J, B)  # or: J × B
end