# Demo - Kuramato-Sivashinsky equation¶

- Authors
Mikael Mortensen (mikaem at math.uio.no)

- Date
Nov 14, 2019

*Summary.* This is a demonstration of how the Python module shenfun can be used to solve the time-dependent,
nonlinear Kuramato-Sivashinsky equation, in a doubly periodic domain. The demo is implemented in
a single Python file KuramatoSivashinsky.py, and it may be run
in parallel using MPI.

## The Kuramato-Sivashinsky equation¶

## Implementation¶

The model equation (1) is implemented in shenfun using Fourier basis functions for both \(x\) and \(y\) directions. We start the solver by implementing necessary functionality from required modules like Numpy, Sympy matplotlib and mpi4py, in addition to shenfun:

```
from sympy import symbols, exp, lambdify
import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
from shenfun import *
```

The size of the problem (in real space) is then specified, before creating
the `TensorProductSpace`

, which is using a tensor product of two Fourier bases as
basis functions. We also
create a `VectorTensorProductSpace`

, since this is required for computing the
gradient of the scalar field `u`

. The gradient is required for the nonlinear
term.

```
# Size of discretization
N = (128, 128)
comm = MPI.COMM_WORLD
K0 = Basis(N[0], 'F', domain=(-30*np.pi, 30*np.pi), dtype='D')
K1 = Basis(N[1], 'F', domain=(-30*np.pi, 30*np.pi), dtype='d')
T = TensorProductSpace(comm, (K0, K1), **{'planner_effort': 'FFTW_MEASURE'})
TV = VectorTensorProductSpace([T, T])
```

Test and trialfunctions are required for assembling the variational forms:

```
u = TrialFunction(T)
v = TestFunction(T)
```

and some arrays are required to hold the solution. We also create an array
`gradu`

, that will be used to compute the gradient in the nonlinear term.
Finally, the wavenumbers are collected in list `K`

. Here one feature is worth
mentioning. The gradient in spectral space can be computed as `1j*K*U_hat`

.
However, since this is an odd derivative, and we are using an even number `N`

for the size of the domain, the highest wavenumber must be set to zero. This is
the purpose of the last keyword argument to `local_wavenumbers`

below.

```
U = Array(T)
U_hat = Function(T)
gradu = Array(TV)
K = np.array(T.local_wavenumbers(True, True, eliminate_highest_freq=True))
```

Note that using this `K`

in computing derivatives has the same effect as
achieved by symmetrizing the Fourier series by replacing the first sum below
with the second when computing odd derivatives.

Here \(\sideset{}{'}\sum\) means that the first and last items in the sum are divided by two. Note that the two sums are equal as they stand (due to aliasing), but only the latter (known as the Fourier interpolant) gives the correct (zero) derivative of the basis with the highest wavenumber.

Sympy is used to generate an initial condition, as stated in Eq (1)

```
# Use sympy to set up initial condition
x, y = symbols("x,y")
ue = exp(-0.01*(x**2+y**2))
ul = lambdify((x, y), ue, 'numpy')
```

Shenfun has a few integrators implemented in the `integrators`

submodule. Two such integrators are the 4th order explicit Runge-Kutta method
`RK4`

, and the exponential 4th order Runge-Kutta method `ETDRK4`

. Both these
integrators need two methods provided by the problem being solved, representing
the linear and nonlinear terms in the problem equation. We define two methods
below, called `LinearRHS`

and `NonlinearRHS`

```
def LinearRHS(self):
# Assemble diagonal bilinear forms
L = -(inner(div(grad(u))+div(grad(div(grad(u)))), v))
return L
def NonlinearRHS(self, U, U_hat, dU):
# Assemble nonlinear term
global gradu
gradu = TV.backward(1j*K*U_hat, gradu)
dU = T.forward(0.5*(gradu[0]*gradu[0]+gradu[1]*gradu[1]), dU)
return -dU
```

The code should, hopefully, be self-explanatory.

All that remains now is to initialize the solution arrays and to setup the
integrator plus some plotting functionality for visualizing the results. Note
that visualization is only nice when running the code in serial. For parallel,
it is recommended to use `HDF5File`

, to store intermediate results to the HDF5
format, for later viewing in, e.g., Paraview.

The solution is initialized as

```
#initialize
X = T.local_mesh(True)
U[:] = ul(*X)
U_hat = T.forward(U, U_hat)
```

And we also create an update function for plotting intermediate results with a cool colormap:

```
# Integrate using an exponential time integrator
plt.figure()
cm = plt.get_cmap('hot')
image = plt.contourf(X[0], X[1], U, 256, cmap=cm)
plt.draw()
plt.pause(1e-6)
count = 0
def update(u, u_hat, t, tstep, **params):
global count
if tstep % params['plot_step'] == 0 and params['plot_step'] > 0:
u = T.backward(u_hat, u)
image.ax.clear()
image.ax.contourf(X[0], X[1], U, 256, cmap=cm)
plt.pause(1e-6)
count += 1
plt.savefig('Kuramato_Sivashinsky_N_{}_{}.png'.format(N[0], count))
```

Now all that remains is to create the integrator and call it

```
if __name__ == '__main__':
par = {'plot_step': 100}
dt = 0.01
end_time = 100
integrator = ETDRK4(T, L=LinearRHS, N=NonlinearRHS, update=update, **par)
#integrator = RK4(T, L=LinearRHS, N=NonlinearRHS, update=update, **par)
integrator.setup(dt)
U_hat = integrator.solve(U, U_hat, dt, (0, end_time))
```