Content provided under a Creative Commons Attribution license, CC-BY 4.0; code under MIT license. (c)2014 Lorena A. Barba, Gilbert Forsyth. Acknowledgement: NSF for support via CAREER award #1149784 and PASI award #1242245. import numpy as np import matplotlib.pyplot as plt %matplotlib inline def sodIC1(): nx = 50 dx = 20./50 CFL = 0.4 c_initial = 374.17 dt = CFL*(dx/c_initial) x = np.linspace(-10,10,50) nt = np.ceil(.01/dt) #initial conditions wl = np.array([1., 0, 100000.]) wr = np.array([0.125, 0, 10000.]) U = np.ones((3,nx)) U[:,:25] = build_U(wl[0],wl[1],wl[2]) U[:,25:] = build_U(wr[0],wr[1],wr[2]) return U, dx, dt, nx, x, int(nt) def sodIC2(): nx = 50 dx = 25./50 CFL = 0.3 c_initial = 374.17 dt = CFL*(dx/c_initial) x = np.linspace(-10,15,50) nt = np.ceil(.01/dt) #initial conditions wl = np.array([1., 0, 100000.]) wr = np.array([0.01, 0, 1000.]) U = np.ones((3,nx)) U[:,:25] = build_U(wl[0],wl[1],wl[2]) U[:,25:] = build_U(wr[0],wr[1],wr[2]) return U, dx, dt, nx, x, int(nt) def build_U(rho, u, p): gamma = 1.4 e = p / ((gamma-1)*rho) e_T = e + u**2/2 U = np.array([[rho],[rho*u],[rho*e_T]]) return U def build_flux(U_in): '''Takes a 3x1 vector U and generates the flux vector for Conserved Euler Eqs''' gamma = 1.4 u1, u2, u3 = U_in[0], U_in[1], U_in[2] F = np.array([u2,\ u2**2/u1+(gamma-1)*(u3-u2**2/(2*u1)),\ u2/u1*(u3+(gamma-1)*(u3-u2**2/(2*u1)))]) return F def build_jacobian(U): rho = U[0,:] u = U[1,:]/U[0,:] gamma = 1.4 p = (gamma-1)*(U[2,:]-.5*U[1,:]**2/U[0,:]) c = np.sqrt(gamma*p/rho) [row, col] = U.shape A = np.zeros((row,row,col)) A[0,:,:] = np.zeros(col), np.ones(col), np.zeros(col) A[1,:,:] = .5*(gamma-3)*u**2, (3-gamma)*u, gamma*np.ones(col)-1 A[2,:,:] = .5*(gamma-2)*u**3-c**2*u/(gamma-1), (3-2*gamma)/2*u**2 + c**2/(gamma-1) return A def laxfriedrichs(U, dx, dt, nx, nt): UN = np.ones((3,nx)) for i in range(nt): UN[:,1:-1] = .5*(U[:,2:]+U[:,:-2])-\ dt/(2*dx)*(build_flux(U[:,2:]) - build_flux(U[:,:-2])) UN[:,0]=UN[:,1] UN[:,-1]=UN[:,-2] U[:]=UN[:] return U U, dx, dt, nx, x, nt = sodIC1() U = laxfriedrichs(U, dx, dt, nx, nt) def decompose_U(U_in): '''Extract pressure, velocity, sound speed, density, entropy and Mach number from U''' gamma = 1.4 vel = U_in[1,:]/U_in[0,:] pres = (gamma - 1)*(U_in[2,:] - .5 * U_in[1,:]**2 / U_in[0,:]) rho = U_in[0,:] c = np.sqrt(gamma * pres / rho) S = pres/rho**gamma return vel, pres, rho, c, S vel, pres, rho, c, S = decompose_U(U) def plot_shock_tube(vel, pres, rho, c, S, x): plt.figure(figsize=(14,7)) plt.subplot(2,3,1) plt.plot(x,vel) plt.title('Velocity') plt.subplot(2,3,2) plt.plot(x,pres) plt.title('Pressure') plt.subplot(2,3,3) plt.plot(x,rho) plt.title('Density') plt.subplot(2,3,4) plt.plot(x,vel/c) plt.title('Mach Number') plt.subplot(2,3,5) plt.plot(x,c) plt.title('Speed of sound') plt.subplot(2,3,6) plt.plot(x,S) plt.title('Entropy') plot_shock_tube(vel, pres, rho, c, S, x) U, dx, dt, nx, x, nt = sodIC2() U = laxfriedrichs(U, dx, dt, nx, nt) vel, pres, rho, c, S = decompose_U(U) plot_shock_tube(vel, pres, rho, c, S, x) def richtmyer(U, dx, dt, nx, nt, damp=0): UN = np.ones((3,nx)) UN_plus = np.ones((3,nx)) UN_minus = np.ones((3,nx)) for i in range(nt): UN_plus[:,:-1] = .5*(U[:,1:]+U[:,:-1]) -\ dt/(2*dx)*(build_flux(U[:,1:]) - build_flux(U[:,:-1])) UN_minus[:,1:] = UN_plus[:,:-1] UN[:,1:-1] = U[:,1:-1] - dt/dx *\ (build_flux(UN_plus[:,1:-1]) - build_flux(UN_minus[:,1:-1])) +\ damp * (U[:,2:] - 2*U[:,1:-1] + U[:,:-2]) UN[:,0] = UN[:,1] UN[:,-1] = UN[:,-2] U[:,:] = UN[:,:] return U U, dx, dt, nx, x, nt = sodIC1() U = richtmyer(U, dx, dt, nx, nt) vel, pres, rho, c, S = decompose_U(U) plot_shock_tube(vel, pres, rho, c, S, x) U, dx, dt, nx, x, nt = sodIC1() U = richtmyer(U, dx, dt, nx, nt, damp=0.1) vel, pres, rho, c, S = decompose_U(U) plot_shock_tube(vel, pres, rho, c, S, x) U, dx, dt, nx, x, nt = sodIC2() U = richtmyer(U, dx, dt, nx, nt) vel, pres, rho, c, S = decompose_U(U) plot_shock_tube(vel, pres, rho, c, S, x) U, dx, dt, nx, x, nt = sodIC2() U = richtmyer(U, dx, dt, nx, nt, damp=0.125) vel, pres, rho, c, S = decompose_U(U) plot_shock_tube(vel, pres, rho, c, S, x) def maccormack(U, dx, dt, nx, nt, damp=0): UN = np.ones((3,nx)) UN_star = np.ones((3,nx)) BC = np.array([U[:,0], U[:,-1]]) for i in range(nt): UN_star[:,1:-1] = U[:,1:-1] - dt/dx *\ (build_flux(U[:,2:]) - build_flux(U[:,1:-1])) +\ damp * (U[:,2:] - 2* U[:,1:-1] + U[:,:-2]) UN_star[:,0] = UN_star[:,1] UN_star[:,-1] = UN_star[:,-2] UN[:,1:-1] = .5 * (U[:,1:-1] + UN_star[:,1:-1]) - dt / (2*dx) *\ (build_flux(UN_star[:,1:-1]) - build_flux(UN_star[:,:-2])) UN[:,0] = UN[:,1] UN[:,-1] = UN[:,-2] U[:,:] = UN[:,:] return U U, dx, dt, nx, x, nt = sodIC1() U = maccormack(U, dx, dt, nx, nt, damp=0) vel, pres, rho, c, S = decompose_U(U) plot_shock_tube(vel, pres, rho, c, S, x) U, dx, dt, nx, x, nt = sodIC1() U = maccormack(U, dx, dt, nx, nt, damp=0.125) vel, pres, rho, c, S = decompose_U(U) plot_shock_tube(vel, pres, rho, c, S, x) U, dx, dt, nx, x, nt = sodIC2() U = maccormack(U, dx, dt, nx, nt, damp=0.125) vel, pres, rho, c, S = decompose_U(U) plot_shock_tube(vel, pres, rho, c, S, x) def sodIC3(): nx = 400 dx = 25./nx CFL = 0.3 c_initial = 374.17 dt = CFL*(dx/c_initial) x = np.linspace(-10,15,nx) nt = np.ceil(.01/dt) #initial conditions wl = np.array([1., 0, 100000.]) wr = np.array([0.01, 0, 1000.]) U = np.ones((3,nx)) U[:,:nx/2] = build_U(wl[0],wl[1],wl[2]) U[:,nx/2:] = build_U(wr[0],wr[1],wr[2]) return U, dx, dt, nx, x, int(nt) U, dx, dt, nx, x, nt = sodIC3() U = maccormack(U, dx, dt, nx, nt, damp=0) vel, pres, rho, c, S = decompose_U(U) plot_shock_tube(vel, pres, rho, c, S, x) U, dx, dt, nx, x, nt = sodIC3() U = maccormack(U, dx, dt, nx, nt, damp=0.05) vel, pres, rho, c, S = decompose_U(U) plot_shock_tube(vel, pres, rho, c, S, x) from IPython.core.display import HTML def css_styling(): styles = open("../styles/custom.css", "r").read() return HTML(styles) css_styling()