# Initialize sympy for symbolic computations import sympy as sy # Mathjax printing with simpy sy.init_printing(use_latex='mathjax') # Declare symbolic variables rho_0, K_0 = sy.symbols('rho_0 K_0', real=True) R = sy.Matrix(2, 2, lambda i,j: 0) Rval = sy.Matrix(2, 1, lambda i,j: 0) # Define symbolic matrix An = sy.Matrix([[0, K_0], [1/rho_0, 0]]) # Calculate eigenvalues and eigenvectors of matrix An # eigvec[i][0] => ith eigenvalue # eigvec[i][2] => ith eigenvector eig_vec = An.eigenvects() # Simplify eigenvalues and eigenvectors and save in symbolic arrays for i in range(2): Rval[i] = sy.simplify(eig_vec[i][0]) for j in range(2): R[j,i] = sy.simplify(eig_vec[i][2][0][j]) # Show eigenvalue and corresponding eigenvectors in a widget from IPython.html.widgets import interact, interactive from IPython.display import display, HTML, Latex # Function to print the eigenvalue and eigenvector def print_eigen(Eigenvalue): display(Latex('Eigenvalue =' + '$' + sy.latex(Rval[Eigenvalue]) + '\hspace{20mm} $' + 'Eigenvector =$' + sy.latex(R[:,Eigenvalue]) + '$')) # Interactive widget to display eigenvalue and eigenvector interact(print_eigen, Eigenvalue= {'one':0,'two':1}); # Declare symbolic variables ZL, ZR, dp, du = sy.symbols('Z_L Z_R Delta_p Delta_u', real=True) R = sy.Matrix(2, 2, lambda i,j: 0) # Define symbolic matrix and vector R = sy.Matrix([[-ZL, ZR], [1, 1]]) dq = sy.Matrix([[dp],[du]]) # Solve system alpha = R.inv()*dq #Simplify and print al = sy.simplify(alpha[0]) ar = sy.simplify(alpha[1]) display(Latex('$\\alpha_L =$' + '$' + sy.latex(al) + '$' )) display(Latex('$\\alpha_R =$' + '$' + sy.latex(ar) + '$' )) # SOLUTION IMPLEMENTATION #import numpy import numpy as np # Set left and right state and deltaq (q = [sig11,sig22,sig12,u,v]) qL = np.array([-1.0, 0.0]) qR = np.array([1.0, 0.0]) dq = np.array([qR[0]-qL[0], qR[1]-qL[1]]) # Set bulk and density (left and right) bulkL = 1; rhoL = 1 bulkR = 4; rhoR = 2 # Define speeds and impedance (left and right) cL = np.sqrt(bulkL/rhoL); ZL = rhoL*cL cR = np.sqrt(bulkR/rhoR); ZR = rhoR*cR # Define the 2 eigenvectors (from columns of Matrix R) r1 = np.array([-ZL, 1]) r2 = np.array([ZR, 1]) # Compute the 2 alphas det = ZL + ZR alL = (-dq[0] + dq[1]*ZR)/det alR = (dq[0] + dq[1]*ZL)/det # Compute middle state qm qm = qL + alL*r1 ## Should be equivalent to #qms = qR - alR*r2 #it is! # Compute waves characteristics for plotting x = np.linspace(-5,5,50) Wc = np.zeros((2,len(x))) Wc[0][:] = -x/cL Wc[1][:] = x/cR #SOLUTION PLOTTING # Required for plotting %matplotlib inline import numpy as np import matplotlib.pyplot as plt from pylab import * from IPython.html.widgets import interact from ipywidgets import StaticInteract, RangeWidget, RadioWidget, DropDownWidget # Plot Riemann solution def plot_Riemann(time): fig = plt.figure(figsize=(15, 5)) # Create subplot for (x,t) plane tmax = 5 ax = fig.add_subplot(1,3,1) ax.set_xlabel('x') ax.set_ylabel('time') ax.axis([min(x),max(x),0,tmax]) # Plot characteristic lines # Acoustic waves in red ax.plot(x,Wc[0][:], '-r', linewidth=2) ax.plot(x,Wc[1][:], '-r', linewidth=2) # Plot time-line in (x,t) plane ax.plot(x, 0*x+time, 'k', linewidth=3) # Create pressure subplot for Riemann solution ax2 = fig.add_subplot(1,3,2) ax2.set_xlabel('x') ax2.set_ylabel('Pressure') ax2.axis([min(x),max(x),-2,2]) # Create Riemann solution vector and plot sol = np.zeros(len(x)) for i in range(len(x)): if x[i] < -cL*time: sol[i] = qL[0] elif x[i] < cR*time: sol[i] = qm[0] else: sol[i] = qR[0] ax2.plot(x,sol, 'k', linewidth = 3) ax2.fill_between(x,-20, sol, facecolor='blue', alpha=0.2) # Create velocity subplot for Riemann solution ax3 = fig.add_subplot(1,3,3) ax3.set_xlabel('x') ax3.set_ylabel('Velocity') ax3.axis([min(x),max(x),-2,2]) # Create Riemann solution vector and plot sol = np.zeros(len(x)) for i in range(len(x)): if x[i] < -cL*time: sol[i] = qL[1] elif x[i] < cR*time: sol[i] = qm[1] else: sol[i] = qR[1] ax3.plot(x,sol, 'k', linewidth = 3) ax3.fill_between(x,-20, sol, facecolor='blue', alpha=0.2) return fig # Create interactive widget to visualize solution #interact(plot_Riemann, time=(0,5,0.1), qvar={'sigma11':0, 'sigma22':1, 'sigma12':2, #'normal velocity u':3, 'transverse velocity v':4}); StaticInteract(plot_Riemann, time=RangeWidget(0,5,0.25)) from IPython.display import HTML HTML('') # NORMAL AND TRANSVERSE RIEMANN SOLVER SOLUTION (cartesian grid) #import numpy import numpy as np # Set left and right state and deltaq (q = [sig11,sig22,sig12,u,v]) qL = np.array([1.0, 0.0, 0.0]) qR = np.array([0.0, 0.0, 0.0]) dq = np.array([qR[0]-qL[0], qR[1]-qL[1],qR[2]-qL[2]]) # Set normal for normal solver nx = 1.0 ny = 0.0 # Set normal for transverse solver nxt = 0.0 nyt = 1.0 # Set grid size (dx2=dx/2) and time step dx2 = 0.01 dy2 = 0.01 dt = 0.01 # Set bulk and density (left and right) bulkL = 1; rhoL = 1 bulkR = 8; rhoR = 2 # Define speeds and impedance (left and right) cL = np.sqrt(bulkL/rhoL); ZL = rhoL*cL cR = np.sqrt(bulkR/rhoR); ZR = rhoR*cR ## NORMAL SOLVER # Define the 2 eigenvectors (from columns of Matrix R) rL = np.array([-ZL, nx, ny]) rR = np.array([ZR, nx, ny]) # Define eigenvalues sL = -cL sR = cR # Compute the 2 alphas det = ZL + ZR alL = (-dq[0] + (nx*dq[1] + ny*dq[2])*ZR)/det alR = (dq[0] + (nx*dq[1] + ny*dq[2])*ZL)/det # Compute wave fluctuations amdq = alL*rL*sL apdq = alR*rR*sR ## TRANSVERSE SOLVER (assuming same material on top and bottom cells) # Define the 2 eigenvectors (from columns of Matrix R) rB = np.array([-ZR, nxt, nyt]) rU = np.array([ZR, nxt, nyt]) # Define eigenvalues sB = -cR sU = cR det = 2.0*ZR beB = (-apdq[0] + (nxt*apdq[1] + nyt*apdq[2])*ZR)/det beU = (apdq[0] + (nxt*apdq[1] + nyt*apdq[2])*ZR)/det # Compute transverse wave fluctuations bmdq = beB*rB*sB bpdq = beU*rU*sU # CREATE INTERACTIVE PLOT # Required for plotting %matplotlib inline import numpy as np import matplotlib.pyplot as plt from matplotlib.path import Path import matplotlib.patches as patches from matplotlib.patches import Ellipse, Polygon from pylab import * from IPython.html.widgets import interact from ipywidgets import StaticInteract, RangeWidget, RadioWidget, DropDownWidget # Plot diagram for transverse solver def plot_trans_diag(dt,transverse_solver_up,transverse_solver_down): x = np.linspace(-2*dx2,2*dx2,50); y1 = 0.0*x + dy2; y2 = 0.0*x - dy2; y = np.linspace(-2*dy2,2*dy2,50); x1 = 0.0*y + dx2; x2 = 0.0*y - dx2; fig = plt.figure() ax = fig.add_subplot(111,aspect='equal') ax.plot(x,y1,'-k',linewidth=2) ax.plot(x,y2,'-k',linewidth=2) ax.plot(x1,y,'-k',linewidth=2) ax.plot(x2,y,'-k',linewidth=2) ax.axis([min(x),max(x),min(y),max(y)]) ax.axis('off') xplus = sR*dt -dx2 xminus = sL*dt -dx2 yplus = sU*dt yminus = sB*dt # Main patches (A+DQ and A-DQ) verts = [-dx2,-dy2], [xplus,-dy2], [xplus,dy2], [-dx2,dy2] poly = patches.Polygon(verts, facecolor='blue', alpha=0.5, linewidth=3) ax.add_patch(poly) verts = [-dx2,-dy2], [xminus,-dy2], [xminus,dy2], [-dx2,dy2] poly = patches.Polygon(verts, facecolor='blue', alpha=0.5, linewidth=3) ax.add_patch(poly) #Transverse patches if (transverse_solver_up=='On'): verts = [-dx2,-dy2], [xplus,-dy2+yplus], [xplus,dy2+yplus], [-dx2,dy2] poly = patches.Polygon(verts, facecolor='yellow', alpha=0.5, linewidth=3) ax.add_patch(poly) if (transverse_solver_down=='On'): verts = [-dx2,-dy2], [xplus,-dy2+yminus], [xplus,dy2+yminus], [-dx2,dy2] poly = patches.Polygon(verts, facecolor='red', alpha=0.5, linewidth=3) ax.add_patch(poly) return fig StaticInteract(plot_trans_diag, dt=RangeWidget(0.0,0.008,0.0004), transverse_solver_up=RadioWidget(['On','Off']), transverse_solver_down=RadioWidget(['On','Off']))