ApproxFunDiffEqSpectral.jl

using ApproxFun, OrdinaryDiffEq, Sundials, BenchmarkTools
using Plots; gr()

S = Fourier()
n = 100 # n = 1000 takes forever
T = ApproxFun.plan_transform(S, n)
Ti = ApproxFun.plan_itransform(S, n)
x = points(S, n)


u₀ = T*cos.(cos.(x-0.1))
D2 = Derivative(S,2)
L = D2[1:n,1:n]

# Simple Heat

function heat(du,u,L,t)
    A_mul_B!(du, L, u)
end
prob = ODEProblem(heat, u₀, (0.0,10.0),L)
@time u = solve(prob, CVODE_BDF(linear_solver=:Diagonal); reltol=1e-8,abstol=1e-8)
plot(x,Ti*u(0.0))
plot!(x,Ti*u(0.5))
plot!(x,Ti*u(2.0))
plot!(x,Ti*u(10.0))

# Simple Reaction

tmp = similar(u₀)
function f(dû,û,tmp,t)
    A_mul_B!(tmp,Ti,û)
    @. tmp = 0.75sqrt(tmp) - tmp
    A_mul_B!(dû,T,tmp)
end

using DiffEqOperators
A = DiffEqArrayOperator(Diagonal(L))

prob = SplitODEProblem(A, f, u₀, (0.0,10.0),tmp)
@time u = solve(prob, Rodas4(autodiff=false))
plot(x,Ti*u(0.0))
plot!(x,Ti*u(0.5))
plot!(x,Ti*u(1.0))
plot!(x,Ti*u(2.0))

@time u = solve(prob, CVODE_BDF())
plot(x,Ti*u(0.0))
plot!(x,Ti*u(0.5))
plot!(x,Ti*u(1.0))
plot!(x,Ti*u(2.0))

@time u = solve(prob, KenCarp4(autodiff=false))
plot(x,Ti*u(0.0))
plot!(x,Ti*u(0.5))
plot!(x,Ti*u(1.0))
plot!(x,Ti*u(2.0))

@time u = solve(prob, LawsonEuler(), dt=0.1)
plot(x,Ti*u(0.0))
plot!(x,Ti*u(0.5))
plot!(x,Ti*u(1.0))
plot!(x,Ti*u(2.0))

@time u = solve(prob, NorsettEuler(), dt=0.1)
plot(x,Ti*u(0.0))
plot!(x,Ti*u(0.5))
plot!(x,Ti*u(1.0))
plot!(x,Ti*u(2.0))

@time u = solve(prob, ETDRK4(), dt=0.4)
plot(x,Ti*u(0.0))
plot!(x,Ti*u(0.5))
plot!(x,Ti*u(1.0))
plot!(x,Ti*u(2.0))

S = Fourier()
n = 100
x = points(S, n)

D2 = Derivative(S,2)[1:n,1:n]
D  = (Derivative(S) → S)[1:n,1:n]
T = ApproxFun.plan_transform!(S, n)
Ti = ApproxFun.plan_itransform(S, n)

u_cache = zeros(100)

p = (D,D2,T,Ti,zeros(100),zeros(100))
û₀ = T*cos.(cos.(x-0.1))
A = 0.1*D2
function burgers(dû,û,p,t)
    D,D2,T,Ti = p
    u = Ti*û
    up = Ti*(D*û)
    dû .= A*û .- T*(u.*up)
end

prob = ODEProblem(burgers, û₀, (0.0,1.0), p)
@time û  = solve(prob, CVODE_BDF(); reltol=1e-5,abstol=1e-5) # 0.6s

plot(x, Ti*û(0.0))
plot!(x, Ti*û(1.0))

Ax = DiffEqArrayOperator(Diagonal(A))
function burgers_f(dû,û,p,t)
    D,D2,T,Ti,u,up = p
    A_mul_B!(u,D,û)
    A_mul_B!(up,Ti,u)
    A_mul_B!(u,Ti,û)
    @. up = u*up
    A_mul_B!(dû,T,up)
    dû .*= -1
end
prob = SplitODEProblem(Ax, burgers_f, û₀, (0.0,1.0), p)

@time û  = solve(prob, CVODE_BDF(); reltol=1e-5,abstol=1e-5) # 0.6s
plot(x, Ti*û(0.0))
plot!(x, Ti*û(1.0))

@time û  = solve(prob, Rodas4(autodiff=false); reltol=1e-5,abstol=1e-5) # 0.6s
plot(x, Ti*û(0.0))
plot!(x, Ti*û(1.0))

@time û = solve(prob, KenCarp4(autodiff=false))
plot(x, Ti*û(0.0))
plot!(x, Ti*û(1.0))

@time û = solve(prob, LawsonEuler(), dt=0.1)
plot(x, Ti*û(0.0))
plot!(x, Ti*û(1.0))

@time û = solve(prob, NorsettEuler(), dt=0.1)
plot(x, Ti*û(0.0))
plot!(x, Ti*û(1.0))

@time û = solve(prob, ETDRK4(), dt=0.1)
plot(x, Ti*û(0.0))
plot!(x, Ti*û(1.0))



# BC Projection Failures
using ApproxFun, OrdinaryDiffEq, Sundials, BenchmarkTools
using Plots; gr()

S = ChebyshevDirichlet{1,1}()
n = 100
D2 = Derivative(S,2)
C = eye(D2)
P = eye(n); P[1:2,1:2] .= 0

function heat(dǔ, ǔ, p, t)
    D2, C, P = p
    dǔ .= D2*ǔ
end

n = 100
T = ApproxFun.plan_transform(S, n)
Ti = ApproxFun.plan_itransform(S, n)
x = points(S, n)

û₀ = T*((1-x.^2).*cos.(x-0.1)+ 0.2)
plot(x, Ti*û₀)

# Project initial condition to boundary
u = Fun(S,û₀)
C=eye(S)[3:end,:]
tmp = [Evaluation(-1);
       Evaluation(1);
       C]\[0.;0.;u]
tmp(1)
û₀2 = tmp.coefficients[1:end-3]
plot(x, Ti*û₀2)

p = (D2[1:n,1:n],C[1:n,1:n],P[1:n,1:n])
prob = ODEProblem(heat, û₀2, (0.0,1.0), p)

# Projection callback
condition(u,t,integrator) = true
function affect!(integrator)
    S = ChebyshevDirichlet{1,1}()
    u = Fun(S,integrator.u)
    C=eye(S)[3:end,:]
    tmp = [Evaluation(-1);
           Evaluation(1);
           C]\[0.;0.;u]
    integrator.u .= @view tmp.coefficients[1:end-3]
end
cb = DiscreteCallback(condition,affect!,save_positions = (false,false))
tstops = collect(0:0.1:1)
@time û  = solve(prob, Rodas4(autodiff=false);
                 tstops = tstops, saveat = tstops,
                 reltol=1e-10,abstol=1e-10, callback = cb)
plot( x, Ti*û[1])
plot!(x, Ti*û[2])
plot!(x, Ti*û[3])
plot!(x, Ti*û[4])
plot!(x, Ti*û[end])

@time û  = solve(prob, Rodas4(autodiff=false);
                 dense = false,
                 reltol=1e-10,abstol=1e-10, callback = cb)
plot( x, Ti*û[1])
plot!(x, Ti*û[2])
plot!(x, Ti*û[3])
plot!(x, Ti*û[4])
plot!(x, Ti*û[end])