%%latex
where f1=1, f2=1/Pe, w(y,0)=0.5e−(y−0.3)2/0.052, w(0,t)=0.
Case 1: Pe=∞ (scalar convection) and Case 2: Pe=103 (scalar convection-diffusion)
import numpy as np import matplotlib.pyplot as plt import scipy from scipy import sparse from scipy.sparse import linalg import random from scipy import interpolate
%matplotlib inline
# Material Parameters Pe = 10.0**3 # Case 1 (No Diffusion) #eps = 0.0 # Case 2 (Convection-Diffusion) eps = 1.0 f1 = 0.7 f2 = eps*1./Pe # Homogeneous Dirichlet boundary conditions (compatible with initial condition) y_l = 0.; y_r = y_l
# spatial interval xL=0.0; Lx=1.5 # time interval T=1.0 # initial time t0=0. # number of time steps n_t = 70 # constant time step size dt =(T-t0)/float(n_t) DT = dt * np.ones(n_t) # number of grid points Nn = 70 # number of cells in grid Ne = Nn-1 # spatial grid size dx = (Lx-xL)/float(Ne) # array of grid points x = np.linspace(xL,xL+Lx,Nn) #where the dirichlet boundaries are # dir_ind = np.array([0,Nn-1]) # out_ind = np.array([],'i')
# Left hand side of 2nd order centrally discretized system main_diag=np.ones((Nn),) diagonals = [(1+ 2*dt*f2/dx**2)*main_diag, (0.5*f1*dt/dx - dt*f2/dx**2)*main_diag, -(0.5*f1*dt/dx + f2*dt/dx**2)*main_diag] offsets = [0, 1, -1] # resetting coefficient matrix for dirichlet boundary diagonals[0][0] = 1.0; diagonals[0][-1]= 1.0 diagonals[1][0] = 0.0 diagonals[2][-2]= 0.0 A = scipy.sparse.diags(diagonals, offsets, shape = (Nn,Nn), format = 'csr') # right hand side for algebraic system b = np.zeros((Nn,),'d')
# initial data x0 = x w0 = 0.5*np.exp(-(x - 0.3)**2/0.05**2)
# Plot initial solution fig = plt.figure() ax = fig.gca() ax.plot(x,w0,'b') plt.xlabel('x') plt.ylabel('w') plt.title('Initial solution')
%%latex The discretized High Fidelity Model is given by
# array of exact time point values tnList = t0*np.ones(n_t+1) for i in np.arange(1,n_t+1): tnList[i] = tnList[i-1]+DT[i-1] t0=tnList[0]; tn = t0; #ii=0; W = np.zeros((w0.shape[0],len(tnList)),'d') X = W.copy() # Lagrangian grid calculations #print X.shape W[:,0] = w0 X[:,0] = x0; for nn,t in enumerate(tnList[1:]): #print "step {0} solving to t={1}".format(nn+1,t) dt=t-tn #allows for variable time step X[:,nn+1] = X[:,0] + (nn+1)*dt*f1 # Lagrangian coordinates b = W[:,nn] #contribution from prev time step b[0] = y_l # left boundary condition b[-1] = y_r # right boundary condition W[:,nn+1] = scipy.sparse.linalg.spsolve(A,b) tn=t
from ipywidgets import interact
def plot_soln(index_t): max_w=np.max(W.flat) min_w=np.min(W.flat) #print 'Max solution is u={0} at x={1}'.format(max_u,np.where(u == max_u)) # print np.where(u.flat == max_u) #print 'Min solution is u={0} at x={1}'.format(min_u,np.where(u.flat == min_u)) # print min_u # print np.where(u.flat == min_u) fig = plt.figure() ax = fig.gca() ax.plot(x,W[:,index_t],'b') plt.xlabel('x') plt.ylabel('w') ax.set_ylim([min_w,max_w]) plt.title('Snapshot at time $t^i=%12.3f$' % (tnList[index_t])) #leg=ax.legend(['$u^i$','$u^{ex}$'],loc=1)
plot_soln(2*n_t/10)
#interact(plot_soln,index_t=(w.shape[1]-1)/2)
Perform SVD for snapshots
S is a diagonal matrix containing the eigenvalues and U is a square matrix containing the eigenvectors of the kernel Kt,
W_svd = W.copy() # The snapshots are centered to have zero mean #ns_tmp = np.ones(n_t+1,'d') #w_mean = np.mean(w,axis=1) #w_svd -= np.outer(w_mean,ns_tmp) U,S,V = np.linalg.svd(W_svd,full_matrices=False) #print np.diag(S)
fig=plt.figure() ax = fig.gca() plt.plot(np.log(S),'b*',label='S') plt.ylabel('singular values of $\mathbf{W}$')
#truncate solution snapshots n_k = 3 U_k = U[:,0:n_k] U_kT= U_k.conj().T print "Truncating to n_k={0}, \sigma={1}".format(n_k,S[n_k]) #Approximation accuracy estimated from sum of squares of singular values acc = np.sum(S[0:n_k+1]**2)/np.sum(S**2) print "Accuracy of reduced representation ={0}".format(acc)
#Rank k approximation Wpod = U_k.dot(U_kT.dot(W))
def plot_pod(index): max_w=np.max(W.flat) min_w=np.min(W.flat) fig = plt.figure() ax = fig.gca() ax.plot(x,W[:,index],'b',x,Wpod[:,index],'r') plt.xlabel('x') plt.ylabel('w') ax.set_ylim([min_w,max_w]) plt.title('solutions at time $t^i=%12.3f$' % (tnList[index])) leg=ax.legend(['$w^i$','$\hat{w}^i$'],loc=1)
plt.plot(x,W[:,0],'k-.',x,Wpod[:,0],'r-.',x,W[:,n_t/3],'k--',x,Wpod[:,n_t/3],'b--',x,W[:,2*n_t/3],'k-',x,Wpod[:,2*n_t/3],'c-',x,W[:,n_t],'k:',x,Wpod[:,n_t],'m:') #interact(plot_pod,index=(0,n_t-1)) plt.title('Comparison of truth solution (black) \n and Eulerian MOR (red, blue, cyan and magenta)')
def plot_pod_error(index): diff = W-Wpod max_err=np.max(diff.flat) min_err=np.min(diff.flat) fig = plt.figure() ax = fig.gca() ax.plot(x,W[:,index]-Wpod[:,index],'k') plt.xlabel('x') plt.ylabel('w') ax.set_ylim([min_err,max_err]) plt.title('POD error at time $t^i=%6.3f$' % (tnList[index])) leg=ax.legend(['$\epsilon_p^i$'],loc=1)
plot_pod_error(n_t/2) #interact(plot_pod_error,index=(0,n_t-1))
Perform SVD for snapshots of Lagrangian coordinates
S^ is a diagonal matrix containing the eigenvalues and U^ is a square matrix containing the eigenvectors of the kernel Kxt,
%%latex
The procedure for the case of a linear 1D scalar convection-diffusion equation involves the following steps -
%%latex
Find optimal reduced bases:
Define:
# Perform SVD for Lagrangian coordinates X_svd = X.copy() # The snapshots are centered to have zero mean #ns_tmp = np.ones(n_t+1,'d') #w_mean = np.mean(w,axis=1) #w_svd -= np.outer(w_mean,ns_tmp) Ux,Sx,Vx = np.linalg.svd(X_svd,full_matrices=False)
fig = plt.figure() ax = fig.gca() plt.plot(np.log(Sx),'b*',label='Sx') plt.ylabel('singular values of $\mathbf{X}$')
# Compute truncated left and right singular matrices for Lagrangian coordinates nx_k = 3 Ux_k = Ux[:,0:nx_k] Sx_full = np.diag(Sx) Vx_k = (Sx_full.dot(Vx))[0:nx_k,:] Ux_kT = Ux_k.conj().T #Full Rank approximation of Lagrangian coordinate snapshots Xpod_full = Ux.dot((Ux.conj().T).dot(X_svd)) #Rank k approximation of Lagrangian coordinate snapshots Xpod_k = Ux_k.dot(Ux_kT.dot(X_svd))
# Preparing input vector for optimization routine S_full = np.diag(S) V_k = (S_full.dot(V))[:n_k,:] Ux_kflat = Ux_k.flatten('F'); Vx_kflat = Vx_k.flatten('F'); U_kflat = U_k.flatten('F'); V_kflat = V_k.flatten('F') print 'dim(Ux)={0}, dim(Vx)={1}, dim(Uw)={2}, dim(Vw)={3}'.format(Nn*nx_k, nx_k*(n_t+1), Nn*n_k, n_k*(n_t+1)) print 'dim(Ux)={0}, dim(Vx)={1}, dim(Uw)={2}, dim(Vw)={3}'.format(Ux_k.shape, Vx_k.shape, U_k.shape, V_k.shape) #print Ux_kflat.shape + Vx_kflat.shape + U_kflat.shape + V_kflat.shape guess_trunc = np.concatenate((Ux_kflat, Vx_kflat, U_kflat, V_kflat),axis = 0) guess_full = np.concatenate((U_kflat, V_kflat),axis = 0) #guess_inp = np.array([Ux_k, Vx_k, U_k, V_k]) print '\nChecking that optimization routine is not indeterminate ....' print 'Flattened full Lagrangian input guess size = {0}'.format(guess_full.shape) print 'Flattened truncated Lagrangian input guess size = {0}'.format(guess_trunc.shape) print 'Flattened residual output size = {0}'.format(len(W.flatten('F')))
''' def residual_err2(tilde,full): W_interp1 = W.copy() W_interp2 = W.copy() if full == 1: # FULL Lagrangian basis Uk = np.reshape(tilde[:Nnn_k],(Nn,n_k),'F') Vk = np.reshape(tilde[Nnn_k:],(n_k,n_t+1),'F') ## For full Lagrangian basis X_tilde = Ux.dot(Sx_full.dot(Vx)) #Xpod_k# elif full == 0: # TRUNCATED Lagrangian basis Uxk = np.reshape(tilde[:Nnnx_k],(Nn,nx_k),'F') Vxk = np.reshape(tilde[Nnnx_k:(Nnnx_k + (n_t+1)nx_k)],(nx_k,n_t+1),'F') Uk = np.reshape(tilde[(Nnnx_k + (n_t+1)nx_k):(Nnnx_k + (n_t+1)nx_k + Nnn_k)],(Nn,n_k),'F') Vk = np.reshape(tilde[(Nnnx_k + (n_t+1)nx_k + Nnn_k):],(n_k,n_t+1),'F') ## For truncated Lagrangian basis X_tilde = Uxk.dot(Vxk)
W_tilde = Uk.dot(Vk) #
#W_err = W - W_tilde
for nn,t in enumerate(tnList[1:]):
f1 = interpolate.interp1d(x0, W[:,nn],fill_value="extrapolate")
f2 = interpolate.interp1d(x0, W_tilde[:,nn],fill_value="extrapolate")
W_interp1[:,nn] = f1(X_tilde[:,nn]) # Solution snapshot in Lagrangian coordinates
W_interp2[:,nn] = f2(X_tilde[:,nn]) # Solution snapshot in Lagrangian coordinates
#w_interp[:,nn] = np.interp(x0, X_tilde[:,nn], w_tilde[:,nn])
return (W_interp1 - W_interp2).flatten('F')
''' from scipy.optimize import leastsq
#full = 1 #lag, mesg = leastsq(residual_err, guess_full, args=(1), xtol=1.49012e-5, ftol=1.49012e-4)
full = 0 lag, mesg = leastsq(residual_err2, guess_trunc, args=(0), xtol=1.49012e-5,ftol=1.49012e-5)
print mesg
''' Ux_kop = np.reshape(lag[:Nnnx_k],(Nn,nx_k),'F') Vx_kop = np.reshape(lag[Nnnx_k:(Nnnx_k + (n_t+1)nx_k)],(nx_k,n_t+1),'F') U_kop = np.reshape(lag[(Nnnx_k + (n_t+1)nx_k):(Nnnx_k + (n_t+1)nx_k + Nnn_k)],(Nn,n_k),'F') #U_kopT = U_kop.conj().T V_kop = np.reshape(lag[(Nnnx_k + (n_t+1)nx_k + Nnn_k):],(n_k,n_t+1),'F')
X_lag = Ux_kop.dot(Vx_kop)
W_temp = U_kop.dot(V_kop)
plt.plot(x,W[:,0],'k-',x,W_temp[:,0],'g-o',x,W[:,n_t/3],'k-',x,W_temp[:,n_t/3],'g-o',x,W[:,2n_t/3],'k-',x,W_temp[:,2n_t/3],'g-o',x,W[:,n_t],'k-',x,W_temp[:,n_t],'g-o') plt.title('Comparison of truth solution (black) \n and Lagrangian MOR (green dots)')
def residual_err(tilde,full): W_interp = W.copy() if full == 1: # FULL Lagrangian basis Uk = np.reshape(tilde[:Nn*n_k],(Nn,n_k),'F') Vk = np.reshape(tilde[Nn*n_k:],(n_k,n_t+1),'F') ## For full Lagrangian basis X_tilde = Ux.dot(Sx_full.dot(Vx)) #Xpod_k# elif full == 0: # TRUNCATED Lagrangian basis Uxk = np.reshape(tilde[:Nn*nx_k],(Nn,nx_k),'F') Vxk = np.reshape(tilde[Nn*nx_k:(Nn*nx_k + (n_t+1)*nx_k)],(nx_k,n_t+1),'F') Uk = np.reshape(tilde[(Nn*nx_k + (n_t+1)*nx_k):(Nn*nx_k + (n_t+1)*nx_k + Nn*n_k)],(Nn,n_k),'F') Vk = np.reshape(tilde[(Nn*nx_k + (n_t+1)*nx_k + Nn*n_k):],(n_k,n_t+1),'F') ## For truncated Lagrangian basis X_tilde = Uxk.dot(Vxk) W_tilde = Uk.dot(Vk) for nn,t in enumerate(tnList[1:]): f = interpolate.interp1d(X_tilde[:,nn], W_tilde[:,nn],fill_value="extrapolate") W_interp[:,nn] = f(x0) #w_interp[:,nn] = np.interp(x0, X_tilde[:,nn], w_tilde[:,nn]) return (W - W_interp).flatten('F')
from scipy.optimize import leastsq # Optimization call with FULL Lagrangian basis #full = 1 #lag, mesg = leastsq(residual_err, guess_full, args=(1), xtol=1.49012e-5, ftol=1.49012e-4) # Optimization call with TRUNCATED Lagrangian basis full = 0 lag, mesg = leastsq(residual_err, guess_trunc, args=(0), xtol=1.49012e-5,ftol=1.49012e-5) print mesg
# Unwrapping output for rank-k full lagrangian basis approximation if full == 1: U_kop = np.reshape(lag[:(Nn*n_k)],(Nn,n_k),'F') #U_kopT = U_kop.conj().T V_kop = np.reshape(lag[(Nn*n_k):],(n_k,n_t+1),'F') # Full Lagrangian snapshot matrix X_lag = Ux.dot(Sx_full.dot(Vx))#X.copy() # Unwrapping output for rank-k optimally truncated lagrangian basis approximation elif full == 0: Ux_kop = np.reshape(lag[:Nn*nx_k],(Nn,nx_k),'F') Vx_kop = np.reshape(lag[Nn*nx_k:(Nn*nx_k + (n_t+1)*nx_k)],(nx_k,n_t+1),'F') U_kop = np.reshape(lag[(Nn*nx_k + (n_t+1)*nx_k):(Nn*nx_k + (n_t+1)*nx_k + Nn*n_k)],(Nn,n_k),'F') #U_kopT = U_kop.conj().T V_kop = np.reshape(lag[(Nn*nx_k + (n_t+1)*nx_k + Nn*n_k):],(n_k,n_t+1),'F') # Optimized truncated Lagrangian coordinate basis X_lag = Ux_kop.dot(Vx_kop) #X_full = Ux.dot(Sx_full.dot(Vx)) #print X_lag.shape
# Small routine to check the residual error function for bugs using the optimized bases #guess_op = np.concatenate((Ux_kop.flatten('F'), Vx_kop.flatten('F'), U_kop.flatten('F'), V_kop.flatten('F')),axis=0) #err = residual_err(guess_op,0) #err = np.reshape(err,(W.shape),'F') #plt.plot(x,err[:,n_t])
#Rank k approximation Wlag = W.copy() # optimized truncated state basis on the Lagrangian grid W_temp = U_kop.dot(V_kop) #w_temp = U_kop.dot(U_kop.conj().T.dot(w)) for nn,t in enumerate(tnList[1:-1]): f = interpolate.interp1d(X_lag[:,nn+1], W_temp[:,nn+1],fill_value="extrapolate") Wlag[:,nn+1] = f(x0) #wlag[:,nn] = np.interp(x0, X_lag[:,nn], w_temp[:,nn])
plt.plot(x,W[:,0],'k-',x,Wlag[:,0],'g-o',x,W[:,n_t/3],'k-',x,Wlag[:,n_t/3],'g-o',x,W[:,2*n_t/3],'k-',x,Wlag[:,2*n_t/3],'g-o',x,W[:,n_t],'k-',x,Wlag[:,n_t],'g-o') plt.title('Comparison of truth solution (black) \n and Lagrangian MOR (green dots)')
plt.plot(x,W[:,0],'k-',x,Wpod[:,0],'r-*',x,Wlag[:,0],'g-o',x,W[:,n_t/2],'k-',x,Wpod[:,n_t/2],'r-*',x,Wlag[:,n_t/2],'g-o',x,W[:,5*n_t/6],'k-',x,Wpod[:,5*n_t/6],'r-*',x,Wlag[:,5*n_t/6],'g-o') #,x,wpod[:,n_t],'k:',x,wlag[:,n_t],'m:' plt.title('Comparison of truth solution (black), Eulerian MOR (red stars) \n and Lagrangian MOR (green dots)')
def plot_lag_error(index): diff = W-Wlag max_err=np.max(diff.flat) min_err=np.min(diff.flat) fig = plt.figure() ax = fig.gca() ax.plot(x,W[:,index]-Wlag[:,index],'k') plt.xlabel('x') plt.ylabel('w') #ax.set_ylim([min_err,max_err]) plt.title('LAG error at time $t^i=%6.3f$' % (tnList[index])) leg=ax.legend(['$\epsilon_p^i$'],loc=1)
plot_lag_error(n_t/4)
# Left hand side of 2nd order centrally discretized system main_diag_online = np.ones((Nn),) diagonals_online = [(1.0 + f2*2.0*dt/(dx**2))*main_diag_online, - f2*dt/(dx**2)*main_diag_online, - f2*dt/(dx**2)*main_diag_online] offsets = [0, 1, -1] #diagonals = [(1+ 2*dt*f2/dx**2)*main_diag_online, (0.5*f1*dt/dx - dt*f2/dx**2)*main_diag_online, -(0.5*f1*dt/dx + f2*dt/dx**2)*main_diag_online] #offsets = [0, 1, -1] # resetting coefficient matrix for dirichlet boundary diagonals_online[0][0] = 1.0; diagonals_online[0][-1]= 1.0 diagonals_online[1][0] = 0.0 diagonals_online[2][-2]= 0.0 A_online = scipy.sparse.diags(diagonals_online, offsets, shape = (Nn,Nn), format = 'csr') # right hand side for algebraic system b_online = np.zeros((Nn,),'d')
# number of time steps nt_online = 100 # constant time step size dt_online =(T-t0)/float(nt_online) DT_online = dt_online * np.ones(nt_online) tnList_online = t0*np.ones(nt_online+1) for i in np.arange(1,nt_online+1): tnList_online[i] = tnList_online[i-1]+DT_online[i-1] #U_Eul_kop = U_kop.copy() #interpolate Lagrangian state basis from Lagrangian to Eulerian coordinates #for nn in np.arange(n_k): # f = interpolate.interp1d(Ux_kop[:,nn], U_kop[:,nn],fill_value="extrapolate") # U_Eul_kop[:,nn] = f(x0) Anew = A.copy() Ahat = U_kop.conj().T.dot(A_online.dot(U_kop)) print Ahat.shape t0 = tnList_online[0]; tn = t0; #ii=0; Z = np.zeros((n_k,len(tnList_online)),'d') Wnew = np.zeros((Nn,len(tnList_online)),'d') print U_kop.conj().T.shape print Z.shape print Wnew.shape Z[:,0] = U_kop.conj().T.dot(w0) Wnew[:,0] = w0 for nn,t in enumerate(tnList_online[1:]): #print "step {0} solving to t={1}".format(nn+1,t) dt=t-tn #allows for variable time step bhat = (Z[:,nn]) #contribution from prev time step bhat[0] = y_l # left boundary condition bhat[-1] = y_r # right boundary condition #print bhat.shape bnew = (Wnew[:,nn]) #contribution from prev time step #bnew[0] = y_l # left boundary condition #bnew[-1] = y_r # right boundary condition Z[:,nn+1] = scipy.sparse.linalg.spsolve(Ahat,bhat) Wnew[:,nn+1] = scipy.sparse.linalg.spsolve(Anew,bnew) tn=t Znew = U_kop.dot(Z) print Znew.shape print X_lag.shape Znew_E = Wnew.copy() #interpolate Lagrangian state basis from Lagrangian to Eulerian coordinates for nn in np.arange(nt_online): f = interpolate.interp1d(X_lag[:,nn], Znew[:,nn],fill_value="extrapolate") Znew_E[:,nn] = f(x0)
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
<ipython-input-39-eed75b7ed6af> in <module>()
53 #interpolate Lagrangian state basis from Lagrangian to Eulerian coordinates
54 for nn in np.arange(nt_online):
---> 55 f = interpolate.interp1d(X_lag[:,nn], Znew[:,nn],fill_value="extrapolate")
56 Znew_E[:,nn] = f(x0)
IndexError: index 71 is out of bounds for axis 1 with size 71
plt.plot(x,Wnew[:,nt_online/5],'k-',x,Znew_E[:,nt_online/5],'g-o') # #,x,Wnew[:,n_t/3],'k-',x,Znew[:,n_t/3],'g-o',x,Wnew[:,2*n_t/3],'k-',x,Znew[:,2*n_t/3],'g-o',x,Wnew[:,n_t],'k-',x,Znew[:,n_t],'g-o' plt.title('Comparison of truth solution (black) \n and Lagrangian MOR (green dots)')
xx = np.arange(0,10); yy = np.array([0, 5, 1, 3, 6, 9, 2, 4, 6, 8])#10*xx+1 # grid point to extrapolate on xxhat = xx[-1] + 10 # Using numpy interpolation with constant extension extrapolation yyhat1 = np.interp(xxhat,xx,yy) # Using scipy interpolation with explicit linear extrapolation F = interpolate.interp1d(xx,yy,fill_value="extrapolate") yyhat2 = F(xxhat) print F(5.6) print 'xhat = {0},yhat1 = {1}'.format(xxhat,yyhat1) print 'xhat = {0},yhat1 = {1}'.format(xxhat,yyhat2) # Plotting original data plt.subplot(2, 1, 1) plt.plot(xx,yy) # Plotting the two types of extrapolated curves plt.subplot(2, 1, 2) plt.plot(np.append(xx,xxhat),np.append(yy,yyhat1),'r',np.append(xx,xxhat),np.append(yy,yyhat2),'g')