Porcupine.jl

Del(resolutions::AbstractVector)
Del(cell::AbstractMatrix)

constructs ∇ operator (derivative, gradient, divergence, curl) using central difference stencil. Because the stencil is of length 3 in each dimension, the result is 2 elements shorter in each dimension than the input. To instead retain the same size, use border=:smooth which pads the input

Example

1d derivative

dx = 0.1
x = 0:dx:.5
y = x .^ 2
d = Del(dx)
@test d(y)≈[ 0.2, 0.4,0.6,0.8]
@test d(y, border=:smooth) ≈ [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]

2d gradient

dy = dx = 0.1
a = [x^2 + y^2 for x in 0:dx:0.5, y in 0:dy:0.5]
∇ = Del([dx, dy])
∇(a)

#=
4×4 Matrix{SVector{2, Float64}}:
 [0.2, 0.2]  [0.2, 0.4]  [0.2, 0.6]  [0.2, 0.8]
 [0.4, 0.2]  [0.4, 0.4]  [0.4, 0.6]  [0.4, 0.8]
 [0.6, 0.2]  [0.6, 0.4]  [0.6, 0.6]  [0.6, 0.8]
 [0.8, 0.2]  [0.8, 0.4]  [0.8, 0.6]  [0.8, 0.8]
=#

2d divergence, curl

curl of 2d Vector field is taken to be vorticity scalar field. In 3d curl produces vorticity vector field.

a = collect.([
    (0, 0) (-1, 0) (0, 0)
    (0, -1) (0, 0) (0, 1)
    (0, 0) (1, 0) (0, 0)
])
∇ = Del([1, 1])
@test ∇ ⋅ a ≈ ∇(a, dot) ≈ [2]'
@test ∇ × a ≈ ∇(a, cross) ≈ [0]'

a = collect.([
    (0, 0) (0, 1) (0, 0)
    (-1, 0) (0, 0) (1, 0)
    (0, 0) (0, -1) (0, 0)
])
∇ = Del([1, 1])
@test ∇ ⋅ a ≈ [0]' 
@test ∇ × a ≈ [-2]' 

Lap(resolutions::AbstractVector)
Lap(cell::AbstractMatrix)

constructs Laplacian operator

# 2d Example
dy = dx = 0.1
a = [x^2 + y^2 for x in 0:dx:0.5, y in 0:dy:0.5]
∇2 = Lap([dx, dy])
∇2(a)

#=
4x4 Matrix{Float64}:
   4.0  4.0  4.0  4.0  
   4.0  4.0  4.0  4.0  
   4.0  4.0  4.0  4.0  
   4.0  4.0  4.0  4.0  
=#