| ''' | |
| Solution of the advection equation by the Crank-Nicolson scheme with | |
| and without added artificial diffusion (dissipation) | |
| by Shelvean Kapita, 2022 | |
| ''' | |
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| from scipy.sparse import diags | |
| from scipy.sparse.linalg import spsolve | |
| def u_0(x): | |
| # val = np.exp(-x**2) | |
| val = np.where(np.abs(x)<0.5, 1 -2*np.abs(x),0) | |
| return val | |
| mu = 0.5 | |
| eps = 1.0 | |
| J = 301 | |
| L = 10 | |
| h = L/(J-1) | |
| k = mu*h | |
| tfinal = 3 | |
| a = -eps/16 | |
| tsteps = int(np.round(tfinal/k)) | |
| x = np.linspace(-2,8,J) | |
| v_0 = u_0(x) | |
| u_ex = u_0(x-tfinal) | |
| d1 = 0.25*mu*np.ones((v_0.size-5,)) | |
| d2 = np.ones((v_0.size-4,)) | |
| d3 = -d1 | |
| diagonals = [d2,d1,d3] | |
| A = diags(diagonals, [0,1,-1]).tocsr() | |
| diagonals1 = [d2,d3,d1] | |
| B = diags(diagonals1,[0,1,-1]).tocsr() | |
| b = B@v_0[2:J-2] | |
| v = np.zeros((v_0.shape)) | |
| u = np.copy(v) | |
| w_0 = np.copy(v_0) | |
| d = a*(w_0[4:J]-4*w_0[3:J-1]+6*w_0[2:J-2]-4*w_0[1:J-3]+w_0[0:J-4]) | |
| f = b+d | |
| for j in range(tsteps): | |
| w = spsolve(A,b) | |
| z = spsolve(A,f) | |
| v[2:J-2] = w | |
| u[2:J-2] = z | |
| v[0],v[1],v[J-1],v[J-2] = 0, 0, 0, 0 | |
| u[0],u[1],u[J-1],u[J-2] = 0, 0, 0, 0 | |
| b = B@w | |
| d = a*(u[4:J]-4*u[3:J-1]+6*u[2:J-2]-4*u[1:J-3]+u[0:J-4]) | |
| f = (B@z)+d | |
| plt.plot(x,v,linewidth='1.5',label='CN',color='blue') | |
| plt.plot(x,u,linewidth='1.5',label='CN+dissipation',color='red') | |
| plt.plot(x,u_ex, linewidth='0.5',color='green',label='Exact') | |
| plt.ylim(-0.5, 2.1) | |
| plt.legend() | |
| plt.show() | |