0

I've tried to make the code below to solve one of the Euler equation for the Gas Dynamics. I understand my code is not perfect, that is why I'm getting an error which I don't understand. The error message is:

" File "eulersys.py", line 50, in <module>

u_data = Evolve_in_One_Timestep(u_data)

File "eulersys.py", line 33, in Evolve_in_One_Timestep

derivative = -0.25*(u[0][k+1]*u[1][k+1] - u[0][k-1]*u[1][k-1])*(dt/dx) + 0.5*sqrt(a)*(u[0][k+1] - 2*u[0][k] + u[0][k-1])*(dt/dx)

IndexError: invalid index to scalar variable."

Any help is welcome.

```
""" """
from numpy import *
from scitools.std import *
import time
import glob,os
for filename in glob.glob('1D_Euler*.png'):
os.remove(filename)
N = 1001; M = 200; dx = 1.0/M; dt = 0.0001644 ; a = 1.0; u = zeros((2,M))
#Create function to set initial values
def Initialize():
K = linspace(0.0,1.0, M+1)
for i in range(M):
if K[i] < 0.5: u[0][i] = 1.0
elif K[i] >= 0.5 and K[i] <= 1 : u[0][i] = 0.125
elif K[i] < 0.5: u[1][i] = 0.0
elif K[i] >= 0.5 and K[i] <= 1 : u[1][i] = 0.0
return u
def Evolve_in_One_Timestep(u):
#create a 2D array of zeros to store the new values of the solution
u_new = zeros((2,M,), float32)
#loop over every gridpoint, except for k=0 and k=M-1
for k in range (1, M-1):
derivative = -0.25*(u[0][k+1]*u[1][k+1] - u[0][k-1]*u[1][k-1])*(dt/dx) + 0.5*sqrt(a)*(u[0][k+1] - 2*u[0][k] + u[0][k-1])*(dt/dx)
u_new[0][k] = u[0][k] + derivative
u_new[0][0] = u_new[0][1]
u_new[0][M-1] = u_new[0][M-2]
return u_new[0]
#Initialize the data array and iterate from t = 0 to the terminal time, T
u_data = Initialize()
counter = 0
for n in range(0,N):
u_data = Evolve_in_One_Timestep(u_data)
t = n*dt
#Plot after every 100 steps
if (n%100==0):
plotlabel= "t = " + str(n*dt) #time labelling
plot(linspace(0.0,1.0,M),u_data,'r0.3',label=plotlabel,legend='t=%4.2f' % t , hardcopy = '1D_Euler%04d.png' % counter )
counter += 1
time.sleep(0.5) # a pause to control movie speed
if __name__=="__main__":
#Make a movie and labels
movie('1D_Euler*.png')
xlabel("Grid points")
ylabel("Values of U for the different time levels")
title("1D Euler equation, nonlinear")
legend()
show()
```