# Initialize sympy for symbolic computations import sympy as sy # Mathjax printing with simpy sy.init_printing(use_latex='mathjax') # Declare symbolic variables nx, ny, rho, lamda, mu = sy.symbols('n_x n_y rho lamda mu', real=True) R = sy.Matrix(5, 5, lambda i,j: 0) Rval = sy.Matrix(5, 1, lambda i,j: 0) # Define symbolic matrix An = -sy.Matrix([[0,0,0,nx*(lamda+2*mu),ny*lamda], [0,0,0,nx*lamda,ny*(lamda+2*mu)], [0,0,0,ny*mu,nx*mu], [nx/rho,0,ny/rho,0,0], [0,ny/rho,nx/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(5): Rval[i] = sy.simplify(eig_vec[i][0].subs(nx**2,1-ny**2).cancel()) for j in range(5): R[j,i] = sy.simplify(eig_vec[i][2][0][j].subs(nx**2,1-ny**2).cancel()) R[j,i] = sy.simplify((ny*R[j,i]).subs(nx**2,1-ny**2).cancel()) R[j,i] = sy.simplify(R[j,i].subs(nx**2+ny**2,1).cancel()) # 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): #print "Eigenvalue, Eigenvector= " #display((Rval[Eigenvalue], R[:,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=(0,4)); #import numpy import numpy as np # Set left and right state and deltaq (q = [sig11,sig22,sig12,u,v]) qL = np.array([0.0, 0.0, 0.0, 4.0, 0.0]) qR = np.array([0.0, 0.0, 0.0, 0.0, 0.0]) dq = np.array([qR[0]-qL[0], qR[1]-qL[1], qR[2]-qL[2], qR[3]-qL[3], qR[4]-qL[4]]) # Set Lame parameters (\mu, \lambda) and density (left and right) muL = 1; muR = 4 lamL = 1; lamR = 4 rhoL = 1; rhoR = 2 # Set normal (unitary) nx = 1.0/np.sqrt(2); ny = 1.0/np.sqrt(2) nx2 = nx*nx ny2 = ny*ny nxy = nx*ny # Define P and S speeds (left and right) bulkL = lamL + 2*muL bulkR = lamR + 2*muR cpL = np.sqrt(bulkL/rhoL) cpR = np.sqrt(bulkR/rhoR) csL = np.sqrt(muL/rhoL) csR = np.sqrt(muR/rhoR) # Define Eigenvalues s = np.array([-cpL, cpR, -csL, csR, 0]) # Define the 4 eigenvectors (from columns of Matrix R) r1 = np.array([lamL + 2*muL*nx2, lamL + 2*muL*ny2, 2*muL*nxy, nx*cpL, ny*cpL]) r2 = np.array([lamR + 2*muR*nx2, lamR + 2*muR*ny2, 2*muR*nxy, -nx*cpR, -ny*cpR]) r3 = np.array([-2*muL*nxy, 2*muL*nxy, muL*(nx2 - ny2), -ny*csL, nx*csL]) r4 = np.array([-2*muR*nxy, 2*muR*nxy, muR*(nx2 - ny2), ny*csR, -nx*csR]) r5 = np.array([ny2, nx2, -nxy, 0, 0]) # Compute the 4 alphas detP = cpR*bulkL + cpL*bulkR detS = csR*muL + csL*muR al1 = (cpR*(dq[0]*nx2 + dq[1]*ny2 + 2*nxy*dq[2]) + bulkR*(nx*dq[3] + ny*dq[4]))/detP al2 = (cpL*(dq[0]*nx2 + dq[1]*ny2 + 2*nxy*dq[2]) - bulkL*(nx*dq[3] + ny*dq[4]))/detP al3 = (csR*(dq[2]*(nx2 - ny2) + nxy*(dq[1] - dq[0])) + muR*(nx*dq[4] - ny*dq[3]))/detS al4 = (csL*(dq[2]*(nx2 - ny2) + nxy*(dq[1] - dq[0])) - muL*(nx*dq[4] - ny*dq[3]))/detS # Compute four middle states qmL, qsL, qsR, qmR --note order of characteristics is s3,s1,s5,s4,s2 # Starting from left side qmL = qL + al1*r1 qsL = qmL + al3*r3 # Starting from right side qmR = qR - al2*r2 qsR = qmR - al4*r4 ## Calculate alpha 5 to test correctness ## Note alpha 5 (al[4]) is not neccesary to calculate when using clawpack, since it doesn't propagate any wave #p1 = cpR*(lamL*(nx2-ny2)*(dq[0]-dq[1]) - 2*muL*(dq[0]*ny2 + dq[1]*nx2) + 4*nxy*dq[2]*(lamL + muL)) #p2 = cpL*(lamR*(nx2-ny2)*(dq[0]-dq[1]) - 2*muR*(dq[0]*ny2 + dq[1]*nx2) + 4*nxy*dq[2]*(lamR + muR)) #p3 = (nx*dq[3] + ny*dq[4])*(2*lamL*muR - 2*lamR*muL) #al5 = -(p1 + p2 + p3)/detP #qsR_test = qsL + al5*r5 # qsR_test must be equal to qsR #print(qsR_test-qsR) # Must be zero or almost zero in all entries # Compute waves characteristics for plotting x = np.linspace(-5,5,50) Wc = np.zeros((len(s),len(x))) Wc[0][:] = x/s[0] Wc[1][:] = x/s[1] Wc[2][:] = x/s[2] Wc[3][:] = x/s[3] # To plot vertical dashed line Wc[4][:] = x/(s[4]+0.000000000001) # 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,variable): vardict = {'sigma 11':0, 'sigma 22':1, 'sigma 12':2, 'normal velocity u':3, 'transverse velocity v':4} qvar = vardict.get(variable) fig = plt.figure(figsize=(10, 4)) # Create subplot for (x,t) plane tmax = 5 ax = fig.add_subplot(1,2,1) ax.set_xlabel('x') ax.set_ylabel('time') ax.axis([min(x),max(x),0,tmax]) # Plot characteristic lines # P-waves in red ax.plot(x,Wc[0][:], '-r', linewidth=2) ax.plot(x,Wc[1][:], '-r', linewidth=2) # S-waves in blue ax.plot(x,Wc[2][:], '-b', linewidth=2) ax.plot(x,Wc[3][:], '-b', linewidth=2) # Non-propagating wave in dashed ax.plot(x,Wc[4][:], '--k', linewidth=2) # Plot time-line in (x,t) plane ax.plot(x, 0*x+time, 'k', linewidth=3) # Create subplot for Riemann solution ax2 = fig.add_subplot(1,2,2) ax2.set_xlabel('x') ax2.set_ylabel('q [' + str(qvar+1) + ']') ax2.axis([min(x),max(x),-max(x),max(x)]) # Create Riemann solution vector and plot sol = np.zeros(len(x)) for i in range(len(x)): if x[i] < -cpL*time: sol[i] = qL[qvar] elif x[i] < -csL*time: sol[i] = qmL[qvar] elif x[i] < 0: sol[i] = qsL[qvar] elif x[i] < csR*time: sol[i] = qsR[qvar] elif x[i] < cpR*time: sol[i] = qmR[qvar] else: sol[i] = qR[qvar] ax2.plot(x,sol, 'k', linewidth = 3) ax2.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), variable=DropDownWidget(['sigma 11','sigma 22','sigma 12', 'normal velocity u','transverse velocity v']))