3d

3d.jl

function literature_example()
    # Problem: Δu + 1000uₓ = f
    # u = 0 on the boundaries
    # f(x, y, z) = exp(xyz) sin(πx) sin(πy) sin(πz)
    # 2nd order central differences (shows serious wiggles)

    # Unknowns per dimension
    N = 50

    # Total number of unknowns
    n = N^3

    # Mesh width
    h = 1.0 / (N + 1)

    # Interior points only
    xs = linspace(0, 1, N + 2)[2 : N + 1]

    # The Laplacian
    Δ = poisson_matrix(Float64, N, 3) ./ -h^2

    # And the dx bit.
    ∂x_1d = spdiagm((fill(-1000.0 / 2h, N - 1), fill(1000.0 / 2h, N - 1)), (-1, 1))
    ∂x = kron(speye(N^2), ∂x_1d)

    # Final matrix and rhs.
    A = Δ + ∂x
    b = reshape([f(x, y, z) for x ∈ xs, y ∈ xs, z ∈ xs], n)

    x, res = bicgstabl(A, b, 2, maxiter = 500, tol = eps())

    # @time xs, reshape(A \ rhs, N, N, N)
end

function poisson_matrix{T}(::Type{T}, n, dims)
    D = second_order_central_diff(T, n);
    A = copy(D);

    for idx = 2 : dims
        A = kron(A, speye(n)) + kron(speye(size(A, 1)), D);
    end

    A
end

# (-1 2 -1) on the diagonal (actuall minus the laplacian)
second_order_central_diff{T}(::Type{T}, dim) = convert(SparseMatrixCSC{T, Int}, SymTridiagonal(fill(2 * one(T), dim), fill(-one(T), dim - 1)))

# For the rhs
f(x, y, z) = exp.(x .* y .* z) .* sin.(π .* x) .* sin.(π .* y) .* sin.(π .* z)