Diffusion Example¶

# dependencies
using LFAToolkit
using LinearAlgebra
using Pkg
Pkg.activate("./")
Pkg.instantiate()
using Plots

# setup
p = 2
dimension = 1

mesh = []
if dimension == 1
   mesh = Mesh1D(1.0)
elseif dimension == 2
   mesh = Mesh2D(1.0, 1.0)
end

# operator
diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)

# compute full operator symbols
numbersteps = 250
maxeigenvalue = 0
θ_min = -π/2
θ_max = 3π/2
θ_step = 2π/(numbersteps-1)
θ_range = θ_min:θ_step:θ_max

# compute and plot smoothing factor
# -- 1D --
if dimension == 1
    # compute
    maxeigenvalues = zeros(numbersteps)

    for i in 1:numbersteps
        θ = [θ_range[i]]
        A = computesymbols(diffusion, θ)
        eigenvalues = [abs(val) for val in eigvals(A)]
        maxeigenvalues[i] = max(eigenvalues...)
        maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])
    end

    # plot
    println("max eigenvalue: ", maxeigenvalue)
    xrange = θ_range/π
    plot(
        xrange,
        xlabel="θ/π",
        xtickfont=font(12, "Courier"),
        maxeigenvalues,
        ytickfont=font(12, "Courier"),
        ylabel="spectral radius",
        linewidth=3,
        legend=:none,
        title="Mass Operator Symbol"
    )
    ylims!(0.0, max(maxeigenvalues...) * 1.1)
# -- 2D --
elseif dimension == 2
    # compute
    maxeigenvalues = zeros(numbersteps, numbersteps)

    for i in 1:numbersteps, j in 1:numbersteps
        θ = [θ_range[i], θ_range[j]]
        A = computesymbols(diffusion, θ)
        eigenvalues = [abs(val) for val in eigvals(A)]
        maxeigenvalues[i, j] = max(eigenvalues...)
        maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])
    end

    # plot
    println("max eigenvalue: ", maxeigenvalue)
    xrange = θ_range/π
    heatmap(
        xrange,
        xlabel="θ/π",
        xtickfont=font(12, "Courier"),
        xrange,
        ylabel="θ/π",
        maxeigenvalues,
        ytickfont=font(12, "Courier"),
        title="Mass Operator Symbol",
        transpose=true,
        aspect_ratio=:equal
    )
    xlims!(θ_min/π, θ_max/π)
    ylims!(θ_min/π, θ_max/π)
end