Views : 15

Comparison of Traditional MOR with Lagrangian MOR for scalar convection and high-Pe scalar convection-diffusion equations

%%latex

$\begin{align} \frac{\partial w}{\partial t} + f_1(y,t,w) \frac{\partial w}{\partial y} &= f_2(y,t,w) \frac{\partial^2 w}{\partial y^2}, \quad \,\forall (y,t)\in [0,1.5]\times[0,1] \\ \frac{d x}{d t} &= f_1(y,t, w) \end{align}$ where $f_1 = 1$ , $f_2 = 1/Pe$ , $w(y,0) = 0.5e^{-(y-0.3)^2/0.05^2}$ , $w(0,t) = 0$ .

Case 1: $Pe = \infty$ (scalar convection) and Case 2: $Pe = 10^3$ (scalar convection-diffusion)

In [1]:

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

In [2]:

%matplotlib inline

In [3]:

# 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

In [4]:

# 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')

In [5]:

# 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')

In [6]:

# initial data
x0 = x
w0 = 0.5*np.exp(-(x - 0.3)**2/0.05**2)

In [7]:

# Plot initial solution
fig = plt.figure()
ax = fig.gca()
ax.plot(x,w0,'b')
plt.xlabel('x')
plt.ylabel('w')
plt.title('Initial solution')

<matplotlib.text.Text at 0x7fd56e07ced0>

Generate Snapshots using High Fidelity Model (HFM)

%%latex The discretized High Fidelity Model is given by $\begin{align} -\left( f_1 \frac{\Delta t}{2 \Delta y} + f_2 \frac{\Delta t}{\Delta y^2}\right) w_{i-1}^{n+1} + \left(1 + 2 f_2 \frac{\Delta t}{\Delta y^2}\right) w_{i}^{n+1} + \left( f_1 \frac{\Delta t}{2 \Delta y} - f_2 \frac{\Delta t}{\Delta y^2}\right) w_{i+1}^{n+1} &= w_i^n \\ x_i^{n+1} &= x_i^n + \Delta t f_1 \end{align}$

In [8]:

# 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

In [9]:

from ipywidgets import interact

In [10]:

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)

In [11]:

plot_soln(2*n_t/10)

In [12]:

#interact(plot_soln,index_t=(w.shape[1]-1)/2)

Perform SVD

Perform SVD for snapshots $\begin{align} \mathbf{W} = \mathbf{USV}^{T} \end{align}$

is a diagonal matrix containing the eigenvalues and is a square matrix containing the eigenvectors of the kernel , $\begin{align} \mathbf{K}_t = \mathbf{WW}^T = \sum_{i=1}^{n_t} \mathbf{w}_i \mathbf{w}_i^T \end{align}$

In [13]:

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)

In [14]:

fig=plt.figure()
ax = fig.gca()
plt.plot(np.log(S),'b*',label='S')
plt.ylabel('singular values of $\mathbf{W}$')

Text(0,0.5,u'singular values of $\\mathbf{W}$')

In [15]:

#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)

Truncating to n_k=3, \sigma=1.80577869635
Accuracy of reduced representation =0.880742079487

In [16]:

#Rank k approximation
Wpod = U_k.dot(U_kT.dot(W))

In [17]:

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)

In [18]:

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)')

Text(0.5,1,u'Comparison of truth solution (black) \n and Eulerian MOR (red, blue, cyan and magenta)')

In [19]:

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)

In [20]:

plot_pod_error(n_t/2)
#interact(plot_pod_error,index=(0,n_t-1))

Perform SVD on Lagrangian coordinate snapshots

Perform SVD for snapshots of Lagrangian coordinates $\begin{align} \mathbf{X} = \mathbf{\hat{U}\hat{S}\hat{V}}^{T} \end{align}$

is a diagonal matrix containing the eigenvalues and is a square matrix containing the eigenvectors of the kernel , $\begin{align} \mathbf{{K}}_{xt} = \mathbf{XX}^T = \sum_{i=1}^{n_t} \mathbf{x}_i \mathbf{x}_i^T \end{align}$

Lagrangian MOR

%%latex

The procedure for the case of a linear 1D scalar convection-diffusion equation involves the following steps - $\begin{itemize} \item The objective is to find the optimal reduced solution snapshot basis (\textbf{in Lagrangian coordinates}), given by $\mathbf{\tilde{U}_w, \tilde{V}_w}$ and the optimal reduced Lagrangian coordinate snapshot basis, given by $\mathbf{{U_x}, {V_x}}$ that minimizes the approximation error for the Eulerian solution snapshots, as given by eq. (\ref{opt}). Here $(\mathbf{\tilde{U}_w\tilde{V}_w})$ is a reduced approximation of the Lagrangian solution snapshots which is mapped to the Eulerian solution snapshots by the linear interpolation operator $\mathcal{P}^x_{\mathbf{U_x,V_x}}$, that uses the reduced Lagrangian coordinate snapshot basis. Finally the approximation error to be minimized is computed by comparing with the Eulerian truth solution $\mathbf{W}$ obtained from the HFM. \item Numerically we first compute the SVD of the Eulerian solution snapshots and the Lagrangian coordinate snapshots to obtain the truncated basis matrices $\mathbf{U_w}$, $\mathbf{V_w}$, $\mathbf{U_x}$ and $\mathbf{V_x}$ as defined in eqs. (\ref{trunc-sol}) and (\ref{trunc-coord}). \item We define a function to compute the pointwise residual error between the truth solution and an approximate Eulerian solution snapshots reconstructed by interpolatiing the reduced Lagrangian solution snapshots on the reduced Lagrangian coordinate basis. This is supplied to the nonlinear least squares optimization routine of numpy as the residual error function input along with a columnwise flattened array of the initial guess for the Lagrangian coordinate snapshots and the Lagrangian solution snapshots, $\mathbf{U_x,V_x},\mathbf{\tilde{U}_w},\mathbf{\tilde{V}_w}$, in that order. The HFM only provides the Eulerian solution snapshots given by $\mathbf{U_w,V_w}$, and these are used instead of $\mathbf{\tilde{U}_w},\mathbf{\tilde{V}_w}$ as the initial guess. \end{itemize}$

%%latex

Find optimal reduced bases: $\begin{align}\label{opt} \min_{\mathbf{U_x, V_x, \tilde{U}_w, \tilde{V}_w}} \, \lvert\lvert \mathbf{W} - \mathcal{P}^x_{\mathbf{U_x,V_x}}(\mathbf{\tilde{U}_w\tilde{V}_w}) \rvert\rvert \end{align}$

Define: $\begin{align} \mathbf{U_w} = [ \mathbf{U} ]_{:,1:n_k} \in \mathbb{R}^{N \times n_k}\, , \mathbf{V_w} = [ \mathbf{SV^T} ]_{1:n_k, :} \in \mathbb{R}^{n_k \times n_t} \label{trunc-sol}\\ \mathbf{U_x} = [ \mathbf{\hat{U}} ]_{:,1:n_{xk}} \in \mathbb{R}^{N \times n_{xk}}\, , \mathbf{V_x} = [ \mathbf{\hat{S}\hat{V}^T} ]_{1:n_{xk}, :} \in \mathbb{R}^{n_{xk} \times n_t} \label{trunc-coord} \end{align}$

In [21]:

# 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)

In [22]:

fig = plt.figure()
ax = fig.gca()
plt.plot(np.log(Sx),'b*',label='Sx')
plt.ylabel('singular values of $\mathbf{X}$')

Text(0,0.5,u'singular values of $\\mathbf{X}$')

In [23]:

# 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))

In [24]:

# 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')))

dim(Ux)=210, dim(Vx)=213, dim(Uw)=210, dim(Vw)=213
dim(Ux)=(70, 3), dim(Vx)=(3, 71), dim(Uw)=(70, 3), dim(Vw)=(3, 71)

Checking that optimization routine is not indeterminate ....
Flattened full Lagrangian input guess size = (423,)
Flattened truncated Lagrangian input guess size = (846,)
Flattened residual output size = 4970

NO USE: Computing a non-optimized approximation by interpolating the reduced Eulerian state basis on the reduced Lagrangian coordinate basis

NO USE: Computing a non-optimized approximation by interpolating the reduced Eulerian state basis on the FULL Lagrangian coordinate basis

Function computing the residual error between the true state snapshot basis and the truncated state snapshot basis interpolated on the truncated Lagrangian basis

''' 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

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_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')

Optimized truncated Lagrangian coordinate basis

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)')

In [25]:

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')

In [26]:

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

2

In [27]:

# 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

In [28]:

# 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])

In [29]:

#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])

In [30]:

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)')

Text(0.5,1,u'Comparison of truth solution (black) \n and Lagrangian MOR (green dots)')

In [31]:

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)')

Text(0.5,1,u'Comparison of truth solution (black), Eulerian MOR (red stars) \n and Lagrangian MOR (green dots)')

In [32]:

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)

In [ ]:

plot_lag_error(n_t/4)

Online evaluation - No idea what's going on here!!!

In [33]:

# 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')

In [39]:

# 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)

(3, 3)
(3, 70)
(3, 101)
(70, 101)
(70, 101)
(70, 71)

---------------------------------------------------------------------------
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

In [ ]:

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)')

A quick demonstration of the difference in "numpy.interp" and "scipy.interpolate.interp1d" routines

In [ ]:

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')