0
$\begingroup$

I'm having trouble with the following implementation of the KS model (see below) found on

Solving numerically the 1D Kuramoto-Sivashinsky Equation using spectral methods

When tN > 300 an overflow occurs on the line:

u_hat2[:,j+1] = (1 / nx) * np.fft.fftshift(np.fft.fft(u[:,j+1]**2))

I tried using a workaround with the convolution theorem with no use, any help would be much apreciated.

import matplotlib.pyplot as plt

from matplotlib.pyplot import cm

nu = 1
L = 100 
nx = 1024

t0 = 0 
tN = 200
dt = 0.05
nt = int((tN - t0) / 0.05)

# wave number mesh
k = np.arange(-nx/2, nx/2, 1)

t = np.linspace(start=t0, stop=tN, num=nt)
x = np.linspace(start=0, stop=L, num=nx)

# solution mesh in real space
u = np.ones((nx, nt))
# solution mesh in Fourier space
u_hat = np.ones((nx, nt), dtype=complex)

u_hat2 = np.ones((nx, nt), dtype=complex)

# initial condition 
u0 = np.cos((2 * np.pi * x) / L) + 0.1 * np.cos((4 * np.pi * x) / L)

# Fourier transform of initial condition
u0_hat = (1 / nx) * np.fft.fftshift(np.fft.fft(u0))

u0_hat2 = (1 / nx) * np.fft.fftshift(np.fft.fft(u0**2))

# set initial condition in real and Fourier mesh
u[:,0] = u0
u_hat[:,0] = u0_hat

u_hat2[:,0] = u0_hat2

# Fourier Transform of the linear operator
FL = (((2 * np.pi) / L) * k) ** 2 - nu * (((2 * np.pi) / L) * k) ** 4
# Fourier Transform of the non-linear operator
FN = - (1 / 2) * ((1j) * ((2 * np.pi) / L) * k)

# resolve EDP in Fourier space
for j in range(0,nt-1):
  uhat_current = u_hat[:,j]
  uhat_current2 = u_hat2[:,j]
  if j == 0:
    uhat_last = u_hat[:,0]
    uhat_last2 = u_hat2[:,0]
  else:
    uhat_last = u_hat[:,j-1]
    uhat_last2 = u_hat2[:,j-1]
  
  # compute solution in Fourier space through a finite difference method
  # Cranck-Nicholson + Adam 
  u_hat[:,j+1] = (1 / (1 - (dt / 2) * FL)) * ( (1 + (dt / 2) * FL) * uhat_current + ( ((3 / 2) * FN) * (uhat_current2) - ((1 / 2) * FN) * (uhat_last2) ) * dt )
  # go back in real space
  u[:,j+1] = np.real(nx * np.fft.ifft(np.fft.ifftshift(u_hat[:,j+1])))
  u_hat2[:,j+1] = (1 / nx) * np.fft.fftshift(np.fft.fft(u[:,j+1]**2))

# plot the result
fig, ax = plt.subplots(figsize=(10,8))

xx, tt = np.meshgrid(x, t)
levels = np.arange(-3, 3, 0.01)
cs = ax.contourf(xx, tt, u.T, cmap=cm.jet)
fig.colorbar(cs)

ax.set_xlabel("x")
ax.set_ylabel("t")
ax.set_title(f"Kuramoto-Sivashinsky: L = {L}, nu = {nu}")```
$\endgroup$
5
  • $\begingroup$ Are the inputs finite? $\endgroup$ Commented Dec 6, 2021 at 13:37
  • $\begingroup$ Yes they are, time goes from t0 to tN, and x goes from 0 to L with periodic boundary conditions $\endgroup$
    – Daniel
    Commented Dec 7, 2021 at 1:05
  • $\begingroup$ I meant the u0 variable in your code, or the u variable in the one line you show above. If u0 is already very (unreasonably) large, then it may not surprise that the FFT overflows. $\endgroup$ Commented Dec 7, 2021 at 17:25
  • $\begingroup$ It is not going to be larger than 1024 elements, u0 is of size nx, which in this case is 1024, but I also tried with only 64 elements and also end up having the overflow $\endgroup$
    – Daniel
    Commented Dec 7, 2021 at 23:07
  • $\begingroup$ No, I was wondering about the size of the elements of that vector, not its length. $\endgroup$ Commented Dec 8, 2021 at 5:09

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.