more examples

examples/fourier1d/ks1.jl

@@ -0,0 +1,86 @@
+#
+using FourierSpaces
+let
+    # add dependencies to env stack
+    pkgpath = dirname(dirname(pathof(FourierSpaces)))
+    tstpath = joinpath(pkgpath, "test")
+    !(tstpath in LOAD_PATH) && push!(LOAD_PATH, tstpath)
+end
+
+using OrdinaryDiffEq, LinearAlgebra
+using Plots, Test
+
+"""
+Kuramoto-Sivashinsky equation
+
+∂ₜu + Δu + Δ²u + 1/2 |∇u|² = 0
+
+x ∈ [0, L)ᵈ (periodic)
+
+TODO: Compute Lyapunov exponent (maybe sensitivities) in 1D/ 2D
+
+https://en.wikipedia.org/wiki/Kuramoto%E2%80%93Sivashinsky_equation
+https://royalsocietypublishing.org/doi/epdf/10.1098/rspa.2014.0932
+"""
+
+N = 256
+T = Float32
+L = T(5pi)
+tspan = (0.0, 10)
+p = nothing
+
+function uIC(x; μ=zero(T), σ=T(0.1))
+    u = @. exp(-1f0/2f0 * ((x-μ)/σ)^2)
+    # reshape(u, :, 1)
+end
+
+""" space discr """
+domain = IntervalDomain(-L, L)
+V = FourierSpace(N; domain = domain)
+Vh = transform(V)
+discr  = Collocation()
+
+(x,) = points(V)
+iftr = transformOp(Vh)
+ftr  = transformOp(V)
+
+û0 = ftr * uIC(x)
+
+function convect!(v, u, p, t)
+    copy!(v, u)
+end
+
+Â = laplaceOp(Vh, discr) # -Δ
+B̂ = biharmonicOp(Vh, discr) # Δ²
+Ĉ = advectionOp((zero(û0),), Vh, discr; vel_update_funcs! =(convect!,)) # uuₓ
+F̂ = SciMLOperators.NullOperator(Vh) # F = 0
+
+L = cache_operator(Â - B̂, û0)
+N = cache_operator(-Ĉ + F̂, û0)
+
+""" time discr """
+tsave = range(tspan...; length=100)
+odealg = Tsit5()
+odealg = SSPRK43()
+prob = SplitODEProblem(L, N, û0, tspan, p)
+
+odecb = begin
+    function affect!(int)
+        if int.iter % 100 == 0
+            println("[$(int.iter)] \t Time $(round(int.t; digits=8))s")
+        end
+    end
+
+    DiscreteCallback((u,t,int) -> true, affect!, save_positions=(false,false))
+end
+@time sol = solve(prob, odealg, saveat=tsave, abstol=1e-4, callback=odecb)
+
+""" analysis """
+pred = iftr.(sol.u, nothing, 0)
+pred = hcat(pred...)
+
+# plot(pred, label=nothing)
+
+anim = animate(pred, V, sol.t, title = "Kuramoto-Sivashinsky 1D")
+gif(anim, joinpath(dirname(@__FILE__), "ks1.gif"), fps=25)
+#

examples/fourier2d/advect.jl

@@ -11,53 +11,77 @@ end
 using OrdinaryDiffEq, LinearAlgebra
 using Test, Plots
 
-nx = 32
-ny = 32
+nx = 128
+ny = 128
 ν = 0e0
 p = nothing
 
 """ space discr """
-V = FourierSpace(nx, ny)
+L = 1.0
+interval = IntervalDomain(-L, L; periodic = true)
+domain = interval × interval
+V = FourierSpace(nx, ny; domain)
 discr = Collocation()
 
 (x,y) = points(V)
 
 """ operators """
 A = -diffusionOp(ν, V, discr)
 
-vx = 1.0; velx = @. x*0 + vx
-vy = 1.0; vely = @. x*0 + vy
+vx = 0.25; velx = @. x*0 + vx
+vy = 0.25; vely = @. x*0 + vy
 C = advectionOp((velx, vely), V, discr)
-F = -C
 
-A = cache_operator(A, x)
-F = cache_operator(F, x)
+odeOp = cache_operator(A - C, x)
+
+function meshplt(
+    x::AbstractMatrix,
+    y::AbstractMatrix,
+    u::AbstractMatrix;
+    a::Real = 45, b::Real = 30, c = :grays, legend = false,
+    kwargs...)
+    plt = plot(x, y, u; c, camera = (a,b), legend, kwargs...)
+    plt = plot!(plt, x', y', u'; c, camera = (a,b), legend, kwargs...)
+end
 
 """ IC """
-uIC(x,y) = @. sin(1x) * sin(1y)
+function uIC(x, y; σ = 0.1, μ = (-0.5, -0.5))
+    r2 = @. (x - μ[1])^2 + (y - μ[2])^2
+    @. exp(-1/2 * r2/(σ^2))
+end
 u0 = uIC(x,y)
 
+u0_re = reshape(u0, nx, ny)
+x_re = reshape(x, nx, ny)
+y_re = reshape(y, nx, ny)
+u0_re = reshape(u0, nx, ny)
+
+# plt = heatmap(u0_re) |> display
+plt = meshplt(x_re, y_re, u0_re)
+display(plt)
+# error("")
+
 """ time discr """
-tspan = (0.0, 10.0)
-tsave = (0, π/4, π/2, 3π/4, 2π,)
+tspan = (0.0, 4.0)
+tsave = LinRange(tspan..., 11)
 odealg = Tsit5()
-prob = SplitODEProblem(A, F, u0, tspan, p)
+prob = ODEProblem(odeOp, u0, tspan, p)
 @time sol = solve(prob, odealg, saveat=tsave, abstol=1e-8, reltol=1e-8)
 
 """ analysis """
-pred = Array(sol)
 
-utrue(x, y, vx, vy, t) = uIC(x .- vx*t, y .- vy*t)
+utrue(x, y, vx, vy, t) = uIC(x .- vx * t, y .- vy * t)
 utr = utrue(x,y,vx,vy,sol.t[1])
 for i=2:length(sol.u)
     ut = utrue(x, y, vx, vy, sol.t[i])
     global utr = hcat(utr, ut)
 end
 
+pred = Array(sol)
 anim = animate(pred, V)
 filename = joinpath(dirname(@__FILE__), "advect" * ".gif")
 gif(anim, filename, fps=5)
 
-err = norm(pred .- utr,Inf)
+err = norm(pred .- utr, Inf)
 @test err < 1e-8
 #