Solution of Euler equation for Gas Dynamics

[Simplicity.]

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()

[woooee]

It is probably one of these
u[0][k+1]
but to find it you will have to put everything on separate lines, usually above the error line so it prints before the message,

derivative = (u[0][k+1]*u[1][k+1] - u[0][k-1]*u[1][k-1])
# and then break it down further, i.e. the last part of the equation
derivative += u[0][k+1]
derivative -= 2*u[0][k]
derivative += u[0][k-1]

and then print 0, k, k+1, k-1, len(u[0]) if that is the line the error message is on. It doesn't matter if the correct result is returned from the above, you just want to find the error.

[Simplicity.]

Thanks a lot woooee for your suggestion.

[raptr_dflo]

Other issues:

In Intialize, the second two elif: blocks will not execute, because you've already handled those ranges in the first two. Instead, move the assignments to indented blocks starting on the following line, and consolidate your u[0][k] and u[1][k] assignments.

Also, since Initialize returns the u array, declare it in there (instead of as a global). Otherwise, don't return anything, and just reference u directly after calling Initialize()

Finally, your error is in the second call of Evolve_in_One_Timestep, since in the first call you return 1-D array u_new[0] and then pass that back in, expecting it to be 2-D. And don't you need to assign values into u[1] as well as u[0]? Or is that a fixed boundary condition?

When you're finished, you need a 1-D array to plot, but instead of trying to return the 1-D array from Evolve..., just pass u_data[0] into plot().

[Simplicity.]

Thank you raptr_dflo, your ideas helped!