src/test/kuramoto.c
Kuramoto–Sivashinsky equation
Duchemin et Eggers, 2014, section 6 propose to use their “Explicit-Implicit-Null” method to solve the Kuramoto–Sivashinsky equation \displaystyle \partial_t u = -u \partial_x u - \partial^2_x u - \partial^4_x u while avoiding the stringent explicit timestep restriction due to the fourth derivative.
In this example we show that this can also be done using the implicit multigrid solver.
#include "grid/multigrid1D.h"
#include "solve.h"
This is the simple explicit discretisation (which is not used).
void solve_explicit (scalar u, double dt)
{
scalar du[];
foreach()
du[] = - u[]*(u[1] - u[-1])/(2.*Delta)
- (u[-1] - 2.*u[] + u[1])/sq(Delta)
- (u[-2] - 4.*u[-1] + 6.*u[] - 4.*u[1] + u[2])/sq(sq(Delta));
foreach()
u[] += dt*du[];
}
int main()
{
We reproduce the same test case as in section 6.2 of Duchemin & Eggers.
init_grid (512);
L0 = 32.*pi;
periodic (right);
scalar u[];
foreach()
u[] = cos(x/16.)*(1. + sin(x/16.));
The timestep is set to 0.1, which is significantly larger than that in Duchemin & Eggers (0.014).
double dt = 1e-1;
// double dt = 1.4e-4;
int i = 0;
TOLERANCE = 1e-6;
for (double t = 0; t <= 150; t += dt, i++) {
if (i % 1 == 0) {
foreach()
fprintf (stdout, "%g %g %g\n", t, x, u[]);
fputs ("\n", stdout);
}
scalar b[];
foreach()
b[] = u[] - dt*u[]*(u[1] - u[-1])/(2.*Delta);
solve (u,
u[] + dt*(u[-1] - 2.*u[] + u[1])/sq(Delta)
+ dt*(u[-2] - 4.*u[-1] + 6.*u[] - 4.*u[1] + u[2])/sq(sq(Delta)),
b);
fprintf (stderr, "%g %d\n", t, solve_stats.i);
// solve_explicit (u, dt);
}
}
The result can be compared to Figure 8 of Duchemin & Eggers.
set term PNG enhanced font ",10"
set output 'sol.png'
set pm3d map
set xlabel 'x'
set ylabel 't'
set xrange [100:0]
set yrange [0:150]
splot 'out' u 2:1:3

Solution of the Kuramoto–Sivashinski equation (script)
References
[duchemin2014] |
Laurent Duchemin and Jens Eggers. The explicit–implicit–null method: Removing the numerical instability of pdes. Journal of Computational Physics, 263:37–52, 2014. [ .pdf ] |