Hello,
I am interested in solving the 1D Kuramoto-Sivashinsky equation for a field u(x,t) with periodic BCs:
\frac{\partial u}{\partial t} = - \frac{\partial^2 u}{\partial x^2} - \frac{\partial^4 u}{\partial x^4} - u \frac{\partial u}{\partial x}, \mbox{ for } x \in [0, L]
I am using the SplitODEProblem structure from OrdinaryDiffEq.jl to solve this semi-linear PDE. I am only seeing fixed-time stepping solvers. Are there adaptive time-stepping available with DiffEq.jl? I am also interested in accelerating this piece of code.
My implementation:
using ApproxFun
using OrdinaryDiffEq
using FFTW
using BandedMatrices
using CairoMakie
#####################
## Setup
struct Kuramoto_Fourier{R,S, U}
L::Float64
Nx::Int64
Δx::Float64
grid::Vector{Float64}
T::R
Tinv::S
D::BandedMatrices.BandedMatrix
D2::BandedMatrices.BandedMatrix
D4::BandedMatrices.BandedMatrix
# Store the diagonal operator -(D^2 + D^4) for the integrating factor
mD24::U
cache::Vector{Float64}
u::Vector{Float64}
end
function Kuramoto_Fourier(L::Float64, Nx::Int64)
Δx = 1/Nx
S = Fourier(0..L)
grid = points(S, Nx)
D2 = Derivative(S,2)[1:Nx,1:Nx]
D = (Derivative(S) → S)[1:Nx,1:Nx]
T = ApproxFun.plan_transform(S, Nx)
Tinv = ApproxFun.plan_itransform(S, Nx)
D4 = Derivative(S,4)[1:Nx,1:Nx]
mD24 = DiffEqArrayOperator(-Diagonal(D4 + D2))
cache = zeros(Nx)
u = zeros(Nx)
return Kuramoto_Fourier{typeof(T), typeof(Tinv), typeof(mD24)}(L, Nx, Δx, grid, T, Tinv, D, D2, D4, mD24, cache, u)
end
# Compute the term -u \partial u /\partial x
function spectral_kuramoto!(dû,û,p::Kuramoto_Fourier,t)
T = p.T
Tinv = p.Tinv
D = p.D
cache = p.cache
u = p.u
mul!(u, D, û)
mul!(cache, Tinv, u)
mul!(u, Tinv, û)
@. cache = cache * u
mul!(u, T, cache)
@. dû = - u
end
#####################
## Example
L = 16*π
Nx = 256
params = Kuramoto_Fourier(L, Nx)
û0 = params.T*map(xi-> 1 + cos(xi/16)*(1 + sin(xi/16)), params.grid)
prob = SplitODEProblem(params.mD24, spectral_kuramoto!, û0, (0.0,300.0), params);
@time sol = solve(prob, ETDRK4(), dt = 0.05)
# Not sure how to use KenCarp58()
# @time sol = solve(prob, KenCarp58(), adaptive =true)
# Snapshots at t = 0 and t = 300
fig = Figure()
ax = Axis(fig[1,1])
lines!(ax, params.grid, params.Tinv*sol(0))
lines!(ax, params.grid, params.Tinv*sol(300))
resize_to_layout!(fig)
fig
# x-t diagram
heatmap(sol.t[1:end], params.grid, hcat(map(ui-> params.Tinv*ui, sol.u[1:end])...)')