explosion of values when simulation Kuramoto Sivashinsky Equation with demo setup

[deephog]

First of all, thank you for creating this amazing library! I have been learning about PDE solvers in the Python environment, and this package really helped me a lot.

One thing I have a concern about is that when I was trying out the KSE solver provided in the demo, I ended up creating a video that has a scale bar that is drastically expanding. I thought it was because of the boundary condition, then I changed the boundary condition to 'periodic' for both axis, the same situation still persisted. Then I set the boundary condition to be constant 0's, the simulation still explodes until it hits 'nan' values.

I copied the code I used, basically copied from the demo with slight changes that don't affect the mechanism. I looked into the source code for clues but I'm not that good at Python so I didn't figure out anything. If you could provide any information that can help me fix this it will be much appreciated!

code:

from pde import DiffusionPDE, PlotTracker, ScalarField, UnitGrid
from pde import PDEBase, ScalarField, UnitGrid
import pde
import pde.visualization.movies as movies
import numpy as np
import time

class KuramotoSivashinskyPDE(PDEBase):
""" Implementation of the normalized Kuramoto–Sivashinsky equation """

def evolution_rate(self, state, t=0):
    """ implement the python version of the evolution equation """
    state_lap = state.laplace(bc=['periodic', 'periodic'])
    state_lap2 = state_lap.laplace(bc=['periodic', 'periodic'])
    state_grad = state.gradient(bc=['periodic', 'periodic'])
    return -state_grad.to_scalar("squared_sum") / 2 - state_lap - state_lap2

def generate_batch(dim=[32, 32], period=6.4, sample_num=64, dt=0.01):
interv = period/sample_num
grid = UnitGrid(dim, periodic=[True, True]) # generate grid
state = ScalarField.random_uniform(grid, seed=int(time.time()%1*100000000)) # generate initial condition
storage = pde.MemoryStorage()
if interv < dt:
raise Exception('sampling interval smaller than dt')

trackers = [
    storage.tracker(interval=interv)
]
eq = KuramotoSivashinskyPDE()  # define the physics
result = eq.solve(state, t_range=period, dt=dt, tracker=trackers)

fields = []
for field in storage:
    fields.append(field.data)
fields = np.asarray(fields)

result.plot()
movies.movie(storage, filename="results/diffusion.mov")
np.save("./results/result.npy", fields)

def main():
generate_batch(dim=[256, 256], period=100, sample_num=100, dt=0.005)
if name == "main":
main()

[david-zwicker]

Could you please post a properly formatted minimal working example, so I can reproduce the problem?

[deephog]

so you do not see a drastically expanding scale bar in the video? The oscillation is getting stronger and stronger almost exponentially, which objects to the energy conservation. I can see this with this exact block of code.

[david-zwicker]

I don't think I see the described behavior with the code listed above. However, I also don't think you should expect the system to be conserved in any way – shouldn't the first term always be negative and thus drive the whole field towards smaller and smaller values?

[gerritwellecke]

I just ran into the same problem and I don't think one can simply see that the equation should always decrease. A recent textbook by Datseris & Parlitz discusses this in its chapter 11.3.2 and mentions that the KS-eq. cannot be solved in real space because of numerical instabilities caused by the fourth spatial derivative. They propose to solve it in Fourier space instead.

[gerritwellecke]

So it turns out that the Kuramoto-Sivashinsky equation exists with different definitions of the squared term. Some use the square of the gradient $|\nabla c|^2$ (the example code does too!) and this equation is not conserved. So the behaviour observed by @deephog is correct. Others use the gradient of the square $\nabla (c^2)$. With periodic boundary conditions this equation seems to be conserved.

This code will run a simulation with the second (i.e. conserved) equation and should yield the desired behaviour according to the OP:

import pde

# make 1D grid
 grid = pde.CartesianGrid([(0, 60)], [128], periodic=True)

# initial condition: random field
state = pde.ScalarField.random_uniform(grid, vmin=-1, vmax=1)

# define Kuramoto-Sivashinksky PDE
# eq = pde.PDE({"u": "-gradient_squared(u) / 2 - laplace(u + laplace(u))"}) <- this was the example
eq = pde.PDE({"u": "-u * d_dx(u) / 2 - laplace(u + laplace(u))"})

# solve the system
storage = pde.MemoryStorage()
result = eq.solve(
    state,
    t_range=50,
    dt=0.001,
    adaptive=True,
    tracker=["progress", storage.tracker(0.1)],
)

# plot the time evolution of the system
pde.plot_kymograph(storage)