Hi all,
I have developed the code below to solve a non-linear scalar hyperbolic conservation law (Burgers' equation). I want to advance to solving a one dimensional system of Euler equations for Gas Dynamics (which is also a hyperbolic partial differential equation) in the similar manner. Can anyone advise on how to do this? Thanks.
""" """
from numpy import *
from scitools.std import *
import time
import glob,os
for filename in glob.glob('Burgers*.png'):
os.remove(filename)
N = 1001
M = 75;
Dx = 1.5/M
Dt = 0.0005
a = 1.0
def Initialise():
U = zeros(M)
K = linspace(0.0,1.5, M+1)
for i in range(M):
if K[i] <= 0.5: U[i] = -0.5
elif K[i] >= 0.5 and K[i] <= 1 : U[i] = 1
elif K[i] >= 1: U[i] = 0
return U
def Evolve_in_One_Timestep (U):
#Create an array to store the new values of the solution
U_new = zeros((M,), float32)
for k in range (1, M-1):
derivative = -0.25*(U[k+1]**2 - U[k-1]**2)*(Dt/Dx) + 0.5*sqrt(a)*(U[k+1] - 2*U[k] + U[k-1])*(Dt/Dx)
U_new[k] = U[k] + derivative
U_new[0] = -0.5
U_new[M-1] = 0
return U_new
#Initialising the data array
U_data = Initialise()
counter = 0
for n in range(0,N):
U_data = Evolve_in_One_Timestep(U_data)
t = n*Dt
if (n%100==0):
plotlabel= "t = " + str(n*Dt) #time labelling
plot(linspace(0.0,0.75,M), U_data,'r0.3',label=plotlabel,legend='t=%4.2f' % t , hardcopy = 'Burgers%04d.png' % counter )
counter += 1
time.sleep(0.5) # a pause to control movie speed
movie('Burgers*.png')
