begin
    import Pkg
    # careful: this is _not_ a reproducible environment
    # activate the local environment
    Pkg.activate(".")
    Pkg.instantiate()
    using PlutoUI, PlutoLinks
    using CairoMakie
end

@revise using Ariadne

Implicit = @ingredients(joinpath(dirname(pathof(Ariadne)), "../examples/implicit.jl"));

using SummationByPartsOperators

function heat_1D!(du, u, (D1m, D1p), t)
    du1 = D1p * u
    mul!(du, D1m, du1)
    return
end

heat_1D! (generic function with 1 method)

# initial conditions
f(x) = sin(x * π)

f (generic function with 1 method)

begin
    xmin = 0.0
    xmax = 1.0

    polydeg = 3
    elements = 40
end

40

begin
    D_local = legendre_derivative_operator(xmin = -1.0, xmax = 1.0, N = polydeg + 1)
    mesh = UniformPeriodicMesh1D(; xmin, xmax, Nx = elements)
end

UniformPeriodicMesh1D{Float64} with 40 cells in (0.0, 1.0)

begin
    D1m = couple_discontinuously(D_local, mesh, Val(:minus))
    D1p = couple_discontinuously(D_local, mesh, Val(:plus))
end

First derivative operator {T=Float64} on 4 Lobatto Legendre nodes in [-1.0, 1.0]
coupled discontinuously (upwind: Val{:plus}()) on UniformPeriodicMesh1D{Float64} with 40 cells in (0.0, 1.0)

x = grid(D1m);

u₀ = f.(x);

lines(x, u₀)

function solve_heat_1D(G!, x, Δt, t_final, initial_condition, p)
    ts = 0.0:Δt:t_final

    u₀ = initial_condition.(x)

    hist = [copy(u₀)]
    callback = (u) -> push!(hist, copy(u))
    Implicit.solve(G!, heat_1D!, u₀, p, Δt, ts; callback)

    return ts, hist
end

solve_heat_1D (generic function with 1 method)

begin
    Δt = 0.01
    t_final = 50.0
end

50.0

ts, hist_E = solve_heat_1D(Implicit.G_Euler!, x, Δt, t_final, f, (D1m, D1p));

_, hist_M = solve_heat_1D(Implicit.G_Midpoint!, x, Δt, t_final, f, (D1m, D1p));

_, hist_T = solve_heat_1D(Implicit.G_Trapezoid!, x, Δt, t_final, f, (D1m, D1p));

@time solve_heat_1D(Implicit.G_Midpoint!, x, Δt, t_final, f, (D1m, D1p));

function plot_timesteps(x, hist, ts, points; title = "")
    fig = Figure()
    ax = Axis(fig[1, 1]; title)
    for p in points
        lines!(ax, x, hist[p], label = "t = $(ts[p])")
    end
    axislegend(ax, position = :cb)
    return fig
end

plot_timesteps (generic function with 1 method)

plot_timesteps(x, hist_E, ts, [1, 2, 3, 4, 5, 10, length(ts)]; title = "Euler Δt=$(Δt)")

plot_timesteps(x, hist_M, ts, [1, 2, 3, 4, 5, 10, length(ts)]; title = "Midpoint Δt=$(Δt)")

plot_timesteps(x, hist_T, ts, [1, 2, 3, 4, 5, 10, length(ts)]; title = "Trapezoid Δt=$(Δt)")

using LinearAlgebra

J_E = Implicit.jacobian(Implicit.G_Euler!, heat_1D!, zero(x), (D1m, D1p), 0.1, 0.0)

160×160 SparseArrays.SparseMatrixCSC{Float64, Int64} with 960 stored entries:
⎡⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⎤
⎢⠀⠈⠿⣧⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠿⣧⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠿⣧⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠘⠻⣦⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠻⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⣀⠀⎥
⎣⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠻⣦⎦

J_M = Implicit.jacobian(Implicit.G_Midpoint!, heat_1D!, zero(x), (D1m, D1p), 0.1, 0.0)

160×160 SparseArrays.SparseMatrixCSC{Float64, Int64} with 960 stored entries:
⎡⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⎤
⎢⠀⠈⠿⣧⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠿⣧⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠿⣧⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠘⠻⣦⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠻⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⣀⠀⎥
⎣⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠻⣦⎦

J_T = Implicit.jacobian(Implicit.G_Trapezoid!, heat_1D!, zero(x), (D1m, D1p), 0.1, 0.0)

160×160 SparseArrays.SparseMatrixCSC{Float64, Int64} with 960 stored entries:
⎡⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⎤
⎢⠀⠈⠿⣧⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠿⣧⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠿⣧⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠘⠻⣦⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠻⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⣀⠀⎥
⎣⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠻⣦⎦

J_M .* 2 + I == J_E

true

J_M == J_T

true

begin
    nnodes = 120
    accuracy_order = 3
    D = upwind_operators(
        periodic_derivative_operator;
        accuracy_order, xmin, xmax, N = nnodes
    )
    x_U = grid(D.minus)
end

0.0:0.008333333333333333:0.9916666666666667

_, hist_EU = solve_heat_1D(Implicit.G_Euler!, x_U, Δt, t_final, f, (D.minus, D.plus));

_, hist_MU = solve_heat_1D(Implicit.G_Midpoint!, x_U, Δt, t_final, f, (D.minus, D.plus));

plot_timesteps(x_U, hist_EU, ts, [1, 2, 3, 4, 5, 10, length(ts)]; title = "Euler-Upwind Δt=$(Δt)")

@bind i PlutoUI.Slider(1:length(ts))

let
    fig = Figure(title = "t = $(ts[i])")
    ax = Axis(fig[1, 1])
    lines!(ax, x, hist_E[i], label = "Euler")
    lines!(ax, x, hist_M[i], label = "Midpoint")
    # lines!(ax, x, hist_T[i], label="Trapezoid")
    axislegend(ax)
    fig
end

lines(x, hist_E[10])

lines(x, hist_M[10])