import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
import funcroe as fc
# Set parameter
N = 1000 # Grid number
cfl = 0.5 # Courant number
xmax = 1
xmin = -1
x = np.linspace(xmin, xmax, N+1) # Position coordinates
dx = (xmax-xmin)/N
n_dis = int(0.5*N)
t_max = 0.14 # Time
gamma = 1.4
# param0 = [0.445, 0.311, 8.928, 0.5, 0.0, 1.4275] # Density Mass Flux Energy
# param0 = [0.345, 0.527, 6.570, 0.5, 0.0, 1.4275]
# param0 = [0.445, 0.311, 8.928, 1.0, 1.994, 7.691]
param0 = [1.0, 0.0, 2.5, 0.25, 0.0, 0.25]
flag = 1
steps = 0
maxsteps = 1e5
current_time = 0
# Set initial value
u_vector = u_0 = fc.euler_generate(N, n_dis, *param0)
u_initial = u_vector
u_t = u_0
# Calculate analytical solution (Riemann sulution)
# initial_3 = [0.445, 0.311, 8.928]
# initial_2 = [0.345, 0.527, 6.570]
# initial_1 = [1.304, 1.994, 7.691]
# initial_0 = [0.5, 0., 1.4275]
# initial_3 = [0.345, 0.527, 6.570]
# initial_2 = [0.345, 0.527, 6.570]
# initial_1 = [1.304, 1.993, 7.688]
# initial_0 = [0.5, 0., 1.4275]
# initial_3 = [0.445, 0.311, 8.928]
# initial_2 = [0.307, 0.580, 5.797]
# initial_1 = [0.943, 1.780, 6.930]
# initial_0 = [1.0, 1.994, 7.691]
initial_3 = [1.0, 0.0, 2.5]
initial_2 = [0.489, 0.386, 1.071]
initial_1 = [0.596, 0.470, 1.104]
initial_0 = [0.25, 0.0, 0.25]
# speed3 = -2.633
# speed2 = -1.636
# speed1 = 1.529
# speed0 = 2.480
# speed3 = -1.636
# speed2 = -1.635
# speed1 = 1.529
# speed0 = 2.479
# speed3 = -2.633
# speed2 = -1.205
# speed1 = 1.889
# speed0 = 3.719
speed3 = -1.183
speed2 = -0.237
speed1 = 0.788
speed0 = 1.358
dis1 = int((speed3*t_max+1.)*N/2.)
dis2 = int((speed2*t_max+1.)*N/2.)
dis3 = int((speed1*t_max+1.)*N/2.)
dis4 = int((speed0*t_max+1.)*N/2.)
rho = np.zeros(N+1)
mass = np.zeros(N+1)
energy = np.zeros(N+1)
# Density
rho[:dis1+1] = initial_3[0]
rho[dis2+1:dis3+1] = initial_2[0]
rho[dis3+1:dis4+1] = initial_1[0]
rho[dis4+1:] = initial_0[0]
rho1 = np.delete(rho,np.linspace(dis1+1,dis2,dis2-dis1).astype('int64'))
# Mass Flux
mass[:dis1+1] = initial_3[1]
mass[dis2+1:dis3+1] = initial_2[1]
mass[dis3+1:dis4+1] = initial_1[1]
mass[dis4+1:] = initial_0[1]
mass1 = np.delete(mass,np.linspace(dis1+1,dis2,dis2-dis1).astype('int64'))
# Energy
energy[:dis1+1] = initial_3[2]
energy[dis2+1:dis3+1] = initial_2[2]
energy[dis3+1:dis4+1] = initial_1[2]
energy[dis4+1:] = initial_0[2]
energy1 = np.delete(energy,np.linspace(dis1+1,dis2,dis2-dis1).astype('int64'))
solution = np.array([rho1, mass1, energy1])
x1 = x[:dis1+1]
x2 = x[dis2+1:dis3+1]
x3 = x[dis3+1:dis4+1]
x4 = x[dis4+1:]
x_solution = np.hstack((x1, x2, x3, x4))
# Calculate numerical solution
while flag:
steps = steps+1
u_0 = fc.consevation_to_physics(u_vector)
c = np.sqrt(gamma*u_0[2, :]/u_0[0, :]) # Speed of sound
dt = cfl*dx/max(abs(u_0[1, :])+c) # Step of time
if steps > maxsteps:
print('warning: maxsteps reach')
break
if current_time+dt >= t_max:
dt = t_max-current_time
flag = 0
current_time = current_time+dt
u_vector = np.column_stack(
(u_vector[:, 0], u_vector[:, 0], u_vector, u_vector[:, -1], u_vector[:, -1]))
dF = fc.fluxDiscre(dx, u_vector) # Calculate dF
u_1 = u_vector[:, 2:-2]-dt*dF
u_1 = np.column_stack((u_1[:, 0], u_1[:, 0], u_1, u_1[:, -1], u_1[:, -1]))
u_2 = 0.75*u_vector[:, 2:-2]+0.25 * \
u_1[:, 2:-2]-0.25*dt*fc.fluxDiscre(dx, u_1)
u_2 = np.column_stack((u_2[:, 0], u_2[:, 0], u_2, u_2[:, -1], u_2[:, -1]))
u_t = 1/3*u_vector[:, 2:-2]+2/3*u_2[:, 2:-2]-2/3*dt*fc.fluxDiscre(dx, u_2)
u_vector = u_t.copy()
# Plot image
fig = plt.figure(figsize=(14,8))
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
ax = plt.subplot(311)
plt.plot(x, u_initial[0, :], color='k', ls=':')
plt.plot(x, u_t[0, :], 'k', ls='--')
plt.plot(x_solution, solution[0, :], 'k')
plt.scatter(x, u_t[0, :], color='', marker='o', edgecolors='k', s=80)
plt.xlim((-1.0,1.0))
plt.ylim((0.2,1.4))
plt.ylabel('Density')
plt.xticks(np.arange(-1.0,1.5,0.5))
xminorLocator = MultipleLocator(0.1)
yminorLocator = MultipleLocator(0.05)
ax.xaxis.set_minor_locator(xminorLocator)
ax.yaxis.set_minor_locator(yminorLocator)
plt.title('Time = '+str(t_max)+' /Roe/')
ax = plt.subplot(312)
plt.plot(x, u_initial[1, :], color='k', ls=':')
plt.plot(x, u_t[1, :], 'k', ls='--')
plt.plot(x_solution, solution[1, :], 'k')
plt.scatter(x, u_t[1, :], color='', marker='o', edgecolors='k', s=80)
plt.xlim((-1.0,1.0))
plt.ylim((-0.5,2.5))
plt.ylabel('Mass Flux')
plt.xticks(np.arange(-1.0,1.5,0.5))
xminorLocator = MultipleLocator(0.1)
yminorLocator = MultipleLocator(0.1)
ax.xaxis.set_minor_locator(xminorLocator)
ax.yaxis.set_minor_locator(yminorLocator)
plt.title('Time = '+str(t_max)+' /Roe/')
ax = plt.subplot(313)
plt.plot(x, u_initial[2, :], color='k', ls=':')
plt.plot(x, u_t[2, :], 'k', ls='--')
plt.plot(x_solution, solution[2, :], 'k')
plt.scatter(x, u_t[2, :], color='', marker='o', edgecolors='k', s=80)
plt.xlim((-1.0,1.0))
plt.ylim((0.0,10))
plt.ylabel('Energy')
plt.xticks(np.arange(-1.0,1.5,0.5))
xminorLocator = MultipleLocator(0.1)
yminorLocator = MultipleLocator(0.5)
ax.xaxis.set_minor_locator(xminorLocator)
ax.yaxis.set_minor_locator(yminorLocator)
plt.title('Time = '+str(t_max)+' /Roe/')
fig.tight_layout()
plt.savefig('./picture/Roe/roe3.eps', format='eps')
plt.show()
评论5