In [1]:
import sys
sys.path.append("..")
from allthethings import PyNetwork, PyPipe_ps
from allthethings import PyBC_opt_dh, PyMystery_BC
import numpy as np
import matplotlib.pyplot as plt
%pylab inline
from matplotlib import rc
In [2]:
#a = 10; these two work with nodf = 16
#fi = "indata/3pipes3.inp" #name of .inp file
#fc = "indata/3pipes3.config" #name of .config file
# a = 120; these two work with ndof = 8(maybe with 16 but for sure not with 15 or 31)
fi = "../indata/3pipes4.inp"
fc = "../indata/3pipes4.config"
mtype = 1 #model used along network edges. 1 for Preissman Slot. 0 for uniform
n0 = PyNetwork(fi, fc, mtype) #a network with Q==0 at boundaries
n1 = PyNetwork(fi, fc, mtype)
print n0
T = n0.T
M = n0.M
s = linspace(0,T,M+1)
L = n0.Ls[2]
q0 = (.9+.3*exp(-((s-T/4.)*10/T)**2))
q1 = 0*q0[0]*np.ones(size(s))
plot(s,q0)
n0.setbVal(2,q0)
n1.setbVal(2,q1)
In [3]:
#n1.showCurrentData()
dt = T/M
V00 = n0.getTotalVolume()
n0.runForwardProblem(dt) #solve up to time T
n1.runForwardProblem(dt)
In [4]:
hi0 = [n0.getHofA(i) for i in range(3)]
hi1 = [n1.getHofA(i) for i in range(3)]
fig,ax = plt.subplots(nrows = n0.Nedges)
x0 = [0,0,100.-n0.Ns[2]]
for k in range(n0.Nedges):
x = np.arange(0,n0.Ls[k], n0.Ls[k]/n0.Ns[k])+x0[k]*np.ones(size(hi1[k]))
ax[k].plot(x,hi0[k],'b')
ax[k].plot(x,hi1[k],'g')
#ax[k].set_ylim([0,ymaxs[i]])
ax[k].set_xlim([0,100])
#ax[k].set_yticks(np.linspace(0,ymaxs[i],6))
In [5]:
pj =0
Nstar = 80
#Nstar = 160
qh0 = n0.qhist(pj)
qh1 = n1.qhist(pj)
p1 = PyPipe_ps(n0.Ns[pj], n0.Ds[pj], n0.Ls[pj], M, n0.a[pj])
In [6]:
def idx_t(i,j,n,N):
return (2*(N+2)*n+(N+2)*i+j)
In [7]:
g = 9.8
hdata0 = np.zeros(M+1)
hdata1 = np.zeros(M+1)
for n in range(M+1):
a0 = qh0[idx_t(0,Nstar,n,n0.Ns[pj])]
a1 = qh1[idx_t(0,Nstar,n,n0.Ns[pj])]
hdata0[n] = p1.Eta(a0,False)/(g*a0)
hdata1[n] = p1.Eta(a1,False)/(g*a1)
In [8]:
figure()
rc('text', usetex=True)
rc('font', family='serif')
plot(s,hdata0,'k-.',linewidth =3)
plot(s,hdata1,'r:', linewidth=2)
xlim(0,T)
legend([r"with $Q_{true}(L,t)$", r"with $Q_i(L,t)$"], loc = 'upper left')
xlabel('$t$ ($s$)')
ylabel('$H^*(t)$ ($m$)')
plt.savefig("../optmysteryhofta120new.eps")
print "H(x*,t) for two scenarios"
#savefig('fig_latex.eps')
In [9]:
xstar =float(Nstar)/n0.Ns[pj]*n0.Ls[pj]
whichnode = 2
ndof =8
modetype = 0
X = q0[0]*ones(ndof)
for i in range(ndof/2):
X[2*i+1] = 0
X = zeros(ndof)
delay = 0
qf = q1
dx = n0.Ls[0]/n0.Ns[0]
a = n0.a[0]
D = n0.Ds[0]
print "xstar=%.2f"%xstar
print "slot width = %.5f m "%(D**2/4*np.pi*9.8/a**2)
print "Optimizing outflow at node %d"%whichnode
In [10]:
opt = PyMystery_BC(fi, fc, ndof, np.array(X), hdata0, modetype,pj, xstar, whichnode,qf, delay)
opt2 = PyMystery_BC(fi, fc, ndof, np.array(X), hdata0, 1,pj, xstar, whichnode,qf, delay)
In [11]:
opt.compute_f()
f0 = opt.f
opt2.compute_f()
f02 = opt2.f
In [12]:
bct = opt.getBCtimeseries(0)
In [13]:
opt.solve()
In [14]:
opt2.solve()
In [15]:
opt.dump()
opt2.dump()
In [16]:
Q1 = opt.getBCtimeseries(0)
Q2 = opt2.getBCtimeseries(0)
ff = opt.f
ff2 = opt2.f
In [17]:
Td = ((n0.Ls[0]-xstar)+n0.Ls[2])/n0.cmax[0]
dM = int(ceil(M*Td/T))
s2 = s[0:M-dM]
semilogy(s,abs(q1-q0),'r',linewidth=2)
semilogy(s, abs(Q1-q0),'b-o')
semilogy(s, abs(Q2-q0),'g--', linewidth = 2)
legend([r"$Q_i(L,t)$, $f$ = %.2e"%f0, r"$Q_f(L,t)$ (Hermite), $f$ = %.2e"%ff,r"$Q_f(L,t)$ (Fourier), $f$ = %.2e"%ff2 ],loc = 'upper left', fontsize =10)
semilogy((T-Td)*ones(10),linspace(1e-7,1000,10),'k:',linewidth=2)
plt.annotate(
'', xy=(T-Td, 6e-6), xycoords='data',
xytext=(T,6e-6), textcoords='data',
arrowprops={'arrowstyle': '<->'})
plt.annotate(
r"$\delta T$", xy=(T-Td+1.5, 1e-5), xycoords='data',
xytext=(5, 0), textcoords='offset points')
xlabel('$t$ ($s$)')
ylabel('$||Q-Q_{true}||$')
ylim(1e-6,1e3)
xlim(0,T)
plt.savefig("../optmysteryerrora120_hfnew.eps")
In [18]:
plot(s, q0[0:],'k-.', linewidth = 3)
plot(s,np.zeros(size(s)),'r:', linewidth =2)
plot(s, Q1[0:],'b', linewidth=2)
plot(s, Q2[0:],'g', linewidth=3)
ylim(-.1,1.3)
xlabel('$t$ [$s$]')
X = opt.x
#plot(linspace(0,50,len(X)/2), X[0::2],'*')
#legend([r"$Q_i(L,t)$", r"$Q_{true}(L,t)$"], loc = 'center right')
legend([r"$Q_{true}(L,t)$",r"$Q_i(L,t)$",r"$Q_f(L,t)$ (Hermite)", "$Q_f(L,t)$ (Fourier)"], loc = 'center left')
plt.annotate(
'', xy=(T-Td, .1), xycoords='data',
xytext=(T,.1 ), textcoords='data',
arrowprops={'arrowstyle': '<->'})
plt.annotate(
r"$\delta T$", xy=(T-Td+1, .15), xycoords='data',
xytext=(5, 0), textcoords='offset points')
ylabel('$Q$ [$m^3/s$]')
xlim(0,T)
plot((T-Td)*ones(10),linspace(-.1,1.3,10),'k:',linewidth=2)
#plt.savefig("../optmysterysoln_a120_hfnew.eps")
Out[18]:
In [19]:
ef = norm(Q1[0:M-dM]-q0[0:M-dM],2)
e0 = norm(q0[0:M-dM])
ef2 = norm(Q2[0:M-dM]-q0[0:M-dM],2)
e02 = norm(q0[0:M-dM])
In [20]:
print "CPU time (s) & %4.2f &%4.2f\\\\ " %(opt.solve_t, opt2.solve_t)
print "actual time (s)& %4.2f& %4.2f \\\\ " % (opt.wsolve_t, opt2.wsolve_t)
print "parallel speedup& %1.1f& %1.1f\\\\ "%(opt.solve_t/opt.wsolve_t, opt2.solve_t/opt2.wsolve_t)
print "$f_i$ & %.4e& %.4e\\\\ "%(f0, f02)
print "$f_f$ & %.4e& %.4e\\\\ " %(ff, ff2)
print "$f_i/f_f$ & %.4e &%.4e \\\\ "%(ff/f0, ff2/f02)
print "$||Q_i - Q_{true}||$ & %.4f& %.4f\\\\" %(e0, e02)
print "$||Q_f - Q_{true}||$ & %.4f &%.4f\\\\" %(ef, ef2)
In [21]:
#plot(s,np.zeros(size(s)),'r:', linewidth =2)
plot(s, q0[0:],'k-.', linewidth = 3)
plot(s,np.zeros(size(s)),'r:', linewidth =2)
ylim(-.1,1.3)
xlabel('$t$ [$s$]')
X = opt.x
legend([r"$Q_{true}(L,t)$", r"$Q_i(L,t)$"], loc = 'center right')
ylabel('$Q$ [$m^3/s$]')
plt.savefig("../optmysteryinitial.eps")