ks implementation using ApproxFun and OrdinaryDiffEq

[JonasIsensee]

Hello there,
I saw that you were interested in an implementation of the Kuramoto-Sivashinsky model
using ApproxFun and OrdinaryDiffEq.
I played around with this for a while but sadly haven't gotten it quite right.
Since I spent a bunch of hours on it, I wanted to share it anyway in case someone else wants to fix it.

Computation is currently slow and solutions decay to 0.
The first can likely be fixed by someone who knows more about these solvers than I do
and the latter seems to be an issue of the wavenumbers k. (resulting in too much diffusion ?)

Note that Fourier() is a real valued sin & cos decomposition.

[ChrisRackauckas]

We'll take a look at it. Thanks for the starter code. There's probably a bit more to do with propagating the linear Diagonal handling.

[MSeeker1340]

Computation is currently slow and solutions decay to 0.
The first can likely be fixed by someone who knows more about these solvers than I do

Yes as @ChrisRackauckas pointed out adding special handling of Diagonal linear operators to the exponential algorithms should help reduce the computation time. Currently Ax is treated as a dense array by the integrators, resulting in inefficient operator caching and mul! operation of the cached operators. This should be easy to fix so let me work on it.

Edit: the discussions above only applies to the caching version of the exponential algorithms (which your test code uses). Krylov versions of the algorithms can already utilize the sparse Diagonal operators and can make computation way faster, especially for large n:

julia> @time sol1 = solve(prob, ETDRK4(); dt=0.25);
  9.364614 seconds (66.41 k allocations: 11.317 GiB, 5.92% gc time)

julia> @time sol2 = solve(prob, ETDRK4(krylov=true, m=20); dt=0.25);
  0.216894 seconds (1.08 M allocations: 97.290 MiB, 6.04% gc time)

and the latter seems to be an issue of the wavenumbers k. (resulting in too much diffusion ?)
Not really sure about this. Running the problem on Tsit5 also makes the answer go to zero so I think this is a problem with the formulation/discretization and not with the exponential integrators.

By the way, implicit integrators currently do not support SplitODEProblem because of interface issues with jacobians, which is a separate issue. However, since the implicit solvers (apart from IMEX) do not make use of the semilinear structure of the right hand side, you can try casting the problem as a regular ODEProblem instead.

[ChrisRackauckas]

By the way, implicit integrators currently do not support SplitODEProblem because of interface issues with jacobians, which is a separate issue. However, since the implicit solvers (apart from IMEX) do not make use of the semilinear structure of the right hand side, you can try casting the problem as a regular ODEProblem instead.

For this we want to test CNAB2. The issue is that it errors since there's a linear operator. We need to put in special handling for the IMEX implicit handling so that way it just solves the linear system when it encounters a linear operator.

[MSeeker1340]

The newest version of ExponentialUtilities.jl (https://github.com/JuliaDiffEq/ExponentialUtilities.jl/releases/tag/v1.1.0) should fix the speed problem for the exponential integrator:

julia> prob, _ = gen_prob(22, 64; dense=false); # original version

julia> prob_dense, _ = gen_prob(22, 64; dense=true); # uses dense `Ax` instead of `Diagonal`

julia> @time sol = solve(prob, ETDRK4(); dt=0.25);
  0.020576 seconds (8.87 k allocations: 935.906 KiB)

julia> @time sol = solve(prob_dense, ETDRK4(); dt=0.25);
  0.293950 seconds (14.82 k allocations: 157.145 MiB, 7.18% gc time)

[ChrisRackauckas]

What timings are we comparing against? I see ~10s in https://github.com/johnfgibson/julia-pde-benchmark/blob/master/1-Kuramoto-Sivashinksy-benchmark.ipynb , so that can't be the right numbers.

When we start really comparing a different integrator though we need to use work-precision diagrams to take into account the difference of error as well. ETDRK4 can probably have an order of magnitude larger dt than CNAB2 and still get the same error for example, so comparisons at the same error will have pretty large differences in dt that should be accounted for. See https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl for how that's done.

But anyways, I think the interesting test is CNAB2. What we want to make sure we can get and show is that a package-based solution to this, using the same algorithm, can match the hand-optimized version. We need a fix to our semilinear handling in our IMEX algorithms (currently our IMEX algorithms are only compatible with nonlinear+nonlinear), but @YingboMa and @MSeeker1340 are working on getting this solved. Since almost all of the time is spent in the linear algebra I don't think this will be a difficult goal to reach, but it will be a milestone for our package ecosystem since the hand-optimized version includes:

• Pseudospectral integration
• Diagonal matrices from the Fourier representation of the derivatives
• An IMEX algorithm
• Specializations on linearity in the implicit handling

and being able to show that we hit all of these features + optimizations is something that will be unique. We're really close!

[MSeeker1340]

What timings are we comparing against? I see ~10s in https://github.com/johnfgibson/julia-pde-benchmark/blob/master/1-Kuramoto-Sivashinksy-benchmark.ipynb , so that can't be the right numbers.

I changed the tspan to match the one used in the testing script (Tmax = 200.0) and for Nx = 10^4 the result is:

julia> prob, _ = gen_prob(10^4/16*pi, 10^4; dense=false);

julia> @time sol = solve(prob, ETDRK4(); dt=1/16);
  6.583057 seconds (1.40 M allocations: 866.530 MiB, 1.65% gc time)

This is about one order of magnitude slower than the benchmark's julia (I assume it uses CNAB2?), but yes as you pointed out with ETDRK4 the time step can be larger, and using Krylov instead of dense caching can also save time.

[ChrisRackauckas]

and using Krylov instead of dense caching can also save time.

Are you sure about that? At this size I assume that dense caching would be faster? The operators are constructed once and then it's an explicit loop with matrix multiplies.

[MSeeker1340]

Ah you're right. My point is correct in regular cases because computing exp(A) for large dense A can be very slow, but for Diagonal A caching is indeed faster.

[JonasIsensee]

I fixed the example.
It now runs nicely showing the typical dynamics.

Beside finding an error in my code,
I encountered this:
JuliaApproximation/ApproxFun.jl#612

which I fixed in my code by defining the derivative matrix myself.