# Adaptive time-stepping for semi-linear ODE (Kuramoto-Sivashinsky)

**URL:** <https://discourse.julialang.org/t/adaptive-time-stepping-for-semi-linear-ode-kuramoto-sivashinsky/108274>\
**Category:** Performance\
**Tags:** question, diffeq, differentialequation\
**Created:** [January 2, 2024, 9:32pm UTC](https://discourse.julialang.org/t/adaptive-time-stepping-for-semi-linear-ode-kuramoto-sivashinsky/108274 "2024-01-02T21:32:13Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [January 2, 2024, 9:32pm UTC](https://discourse.julialang.org/t/adaptive-time-stepping-for-semi-linear-ode-kuramoto-sivashinsky/108274/1 "2024-01-02T21:32:13Z")

</div>

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:

```julia
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])...)')

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 3, 2024, 12:19pm UTC](https://discourse.julialang.org/t/adaptive-time-stepping-for-semi-linear-ode-kuramoto-sivashinsky/108274/2 "2024-01-03T12:19:49Z")

</div>

> [@mleprovost](#):
>
> I am using the `SplitODEProblem` structure from [OrdinaryDiffEq.jl](https://juliahub.com/ui/Packages/OrdinaryDiffEq) to solve this semi-linear PDE. I am only seeing fixed-time stepping solvers.

All of the additive RK methods are adaptive. See the documentation:

> **[Split ODE Solvers · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/split_ode_solve/#OrdinaryDiffEq.jl)**
>
> Documentation for DifferentialEquations.jl.

> [@mleprovost](#):
>
> ```julia
> mD24 = DiffEqArrayOperator(-Diagonal(D4 + D2))
> 
> ```

This was removed years ago. Where are you finding this?

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [January 3, 2024, 2:16pm UTC](https://discourse.julialang.org/t/adaptive-time-stepping-for-semi-linear-ode-kuramoto-sivashinsky/108274/3 "2024-01-03T14:16:39Z")

</div>

Thank you.

The `mD24 = DiffEqArrayOperator(-Diagonal(D4 + D2))` piece comes from some early pandemic code. What would be the modern way of writing the linear operator?

I tried to use the solver `KenCarp3() but I get the warning message

```julia
sol = solve(prob, KenCarp3(), adaptive =true)

```

```julia
┌ Warning: Split ODE problem do not work with factorization linear solvers. Bug detailed in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1643. Defaulting to linsolve=KrylovJL()
└ @ OrdinaryDiffEq /state/partition1/llgrid/pkg/julia/julia-1.9.2/local/share/julia/packages/OrdinaryDiffEq/SUoOD/src/alg_utils.jl:248

```

How to address this?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 3, 2024, 2:34pm UTC](https://discourse.julialang.org/t/adaptive-time-stepping-for-semi-linear-ode-kuramoto-sivashinsky/108274/4 "2024-01-03T14:34:33Z")

</div>

Change to GMRES. I’ve been wanting to handle that for a bit but haven’t gotten to it.