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()
3
Contributors
4
Replies
5
Views
6 Years
Discussion Span
Last Post by Simplicity.
0

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.

Edited by woooee: n/a

0

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

This topic has been dead for over six months. Start a new discussion instead.
Have something to contribute to this discussion? Please be thoughtful, detailed and courteous, and be sure to adhere to our posting rules.