In [5]:
import sys
sys.path.append("..")
from allthethings import PyNetwork, PyPipe_ps
from allthethings import PyBC_opt_dh
import numpy as np
import matplotlib.pyplot as plt
from writeit import rewritePipes
%pylab inline
In [21]:
Hbar = []
fi_old ="../indata/loop.inp"
fc_old = "../indata/loop.config"
mtype =1
Np=8
Mrs = [0.0]*Np
Ds = [1.]*Np
jt = [1,3,2,2,2,2,3,1]
bt = [1,1,1,1,1,1,1,1]
bv = [0,0,0,0,0,0,0,0]
r = [0,0,0,0,0,0,0,1]
h0s =[1,1,1,1,1,1,1,1]
q0s =[1,1,1,1,1,1,1,1]
Ls = [100.,100,100,100,100.,100.,100.,100.]
Ns = [int(L) for L in Ls]
T = 50
M = 6000
Mi =10
a = 100
for k in range(5):
K = 2*k
elevs =np.linspace(0,K,len(jt))
fn = "../indata/newloop"
oldinp = "../indata/loop.inp"
(fi, fc) = rewritePipes(fn,oldinp, Ns, Ls, Mrs, Ds, jt, bt, bv, r, h0s, q0s, T, M, a,elevs)
n0 = PyNetwork(fi, fc, mtype)
M = n0.M
dt = n0.T/n0.M
n0.runForwardProblem(dt) #solve up to time T
Hbar.append([n0.getAveGradH(i) for i in range(n0.M+1)])
In [22]:
from matplotlib import cm
import matplotlib.colors as colors
cmap = cm.get_cmap('ocean')
cnorm =colors.Normalize(vmin=0, vmax = len(Hbar)+1)
scalarMap = cm.ScalarMappable(norm =cnorm, cmap= cmap)
t = np.linspace(0,n0.T,n0.M+1)
for kk in range(len(Hbar)):
cval = scalarMap.to_rgba(kk)
plot(t,Hbar[kk],color = cval, label="%d"%kk)
#print "%f %f" %(max(Hbar[kk]),mean(Hbar[kk]))
#legend(["max %f , mean %f" %(max(Hbar),mean(Hbar))], loc ='upper left')
legend()
gca().set_ylabel(r'$\langle dH/dx\rangle$')
Out[22]:
In [8]:
print elevs
In [9]:
hi = [n0.getHofA(i) for i in range(n0.Nedges)]
fig,ax = plt.subplots(nrows = n0.Nedges, figsize=(15,15))
x00 = [0,0,0]
for k in range(n0.Nedges):
x = np.arange(0,n0.Ls[k], n0.Ls[k]/n0.Ns[k])
ax[k].plot(x,hi[k],'b')
#ax[k].set_ylim((0,.5))
#ax.set_xticklabels([])
#ax[k].set_yticklabels([])
legend(['Q=0', 'reflect'])
print n0.cmax
print dt/(n0.Ls[0]/n0.Ns[0])
print dt
dx = n0.Ls[0]/n0.Ns[0]
In [10]:
print "CFL = %f " %(dt/dx*n0.a[0])
Hbar = [n0.getAveGradH(i) for i in range(n0.M+1)]
t = np.linspace(0,n0.T,n0.M+1)
kk = 1
plot(t,Hbar,'m')
legend(["max %f , mean %f" %(max(Hbar),mean(Hbar))], loc ='upper left')
#savefig("../../gradH_r66.pdf", format='pdf')
title("L0/L2 = %.f/%.f and M = %d, dx = %f"%(n0.Ls[0],n0.Ls[2], n0.M, n0.Ls[0]/n0.Ns[0]))
print "space-average grad(H)"
print " Q = 0 reflecting"
print "max %f " %(max(Hbar)/sqrt(n0.M))
print "mean %f " %(mean(Hbar))
dx = n0.Ls[0]/n0.Ns[0]
print n0.cmax
print n0.a
In [11]:
def idx_t(i,j,n,N):
return (2*(N+2)*n+(N+2)*i+j)
In [12]:
qh = [n0.qhist(i) for i in range(n0.Nedges)]
j = 20#n0.Ns[0]
K = 2
p1 = PyPipe_ps(n0.Ns[0], n0.Ds[K], n0.Ls[K], n0.M, n0.a[K])
qtest = [p1.HofA(qh[K][idx_t(0,j,n,n0.Ns[K])],False) for n in range(n0.M+1)]
plot(t,qtest)
Af = np.pi/4.
print (Af*9.8)/(n0.a[0]**2)
#plot(t,Q00)
In [ ]:
In [13]:
wn=0;
opt1 = PyBC_opt_dh(fi, fc, ndof, x0, wn, Vin, md)
Q0orig = opt1.getBCtimeseries()
opt1.compute_f()
In [ ]:
plot(Q0orig)
Vin2 = np.trapz(Q0orig, x=t)
print Vin2-Vin
print Vin
print "Vin2 = %f"%Vin2
print Q0orig[0]
Dt = T/((ndof-1)/2.)
print Dt
print T/7.*17.5-Vin
print Q0orig[0]
print Vin/T
In [11]:
f0 = opt1.f
print "f is %f" %f0
print "T is %f" %opt1.T
print "Using %s modes" %opt1.modetype
In [12]:
opt1.solve()
In [13]:
print opt1.f
#opt2 = PyBC_opt_dh(fi, fc, 15, np.array(xgood2), wn, Vin)
#opt2.compute_f()
#xgood = opt1.x
In [14]:
ff = opt1.f
#print ff/f0
#print xgood
print opt1.x
#print 2*Vin/n0.T
In [15]:
#make a place to put the data
import os
lucky = 1e6 #you feeling lucky today pal? 1/lucky are your odds of overwriting old data
kk = np.int(lucky*np.random.rand())
where = "../results_ndof%d_wavespeed_%03.f_improvement_%f/"%(ndof,n0.a[0],ff/f0)
if not os.path.exists(where):
os.makedirs(where)
In [16]:
rc('text', usetex=True)
rc('font', family='serif')
ff = opt1.f
Q01 = opt1.getBCtimeseries()
#Q02 = opt2.getBCtimeseries()
t = np.linspace(0,opt1.T, opt1.M+1)
fig,ax = plt.subplots(nrows = 1)
ax.plot(t,Q0orig,'k:', linewidth = 2)#,t,np.zeros(opt1.M+1),'r')
ax.plot(t,Q01,'g', linewidth = 2)
ax.set_ylabel('Q1(t)')
#ax[0].plot(t,Q1,'b')#,t,np.zeros(opt1.M+1),'r')
#ax[0].plot(t,Q12,'g')
#ax[0].set_ylabel('Q1(t)')
#ax[1].plot(t,Q2,'b')#,t,np.zeros(opt1.M+1),'r')
#ax[1].plot(t,Q22, 'g')
s0 = "%.3e"%f0
sf = "%.3e"%ff
w0 = s0.find('e')
wf = sf.find('e')
p0 = int(s0[w0+2:])
pf = int(sf[wf+2:])
legend([r'original, $f = %s\times10^%1.d$'%(s0[0:w0],p0), r'optimized,$f = %s\times10^%1.d$'%(sf[0:wf],pf)], loc = 'upper right')
ax.set_xlabel('t (s)')
#ax[1].set_xlabel('t')
#ax[1].set_ylabel('Q2(t)')
savefig(where+"opt_dh_Qblahnod%d.eps"%ndof, format='eps')
#print max(abs(Q2[i]) for i in range(len(Q2)))
print opt1.x
#print opt2.x
#print "max discrepancy is %e"%max([Q1[i]-Q12[i] for i in range(len(Q1))])
In [17]:
n2 = PyNetwork(fi, fc, mtype)
n3 = PyNetwork(fi,fc, mtype)
n2.setbVal(0,Q0orig)
n3.setbVal(0,Q01)
n2.runForwardProblem(dt)
n3.runForwardProblem(dt)
In [18]:
xf = opt1.x
print xf
In [19]:
Hbar2 = [n2.getAveGradH(i) for i in range(n0.M+1)]
Hbar3 = [n3.getAveGradH(i) for i in range(n0.M+1)]
plot(t,Hbar2, 'k', t,Hbar3,'g', linewidth=2)
ax = gca()
legend(['original', 'optimized'], loc = 'upper left')
ax.set_xlabel('t (s)')
ax.set_ylabel(r"$\langle dH \rangle$")
savefig(where+"opt_dh_dhblahnodf%d.eps"%ndof, format='eps')
In [20]:
with open(where+"more_info.txt",'w') as fout:
fout.write("data (Q0, Qf) pickled to Q0.p and Qf.p\n\n")
fout.write("to generate whole history, load them and run:\n\
n0 = PyNetwork(fi,fc, 1)\n\
n0.setbVal(0,Q)\n\
n0.runForwardProblem(dt)\n\
for Q = Q00 or Q = Qf\n\n")
fout.write("***************************************************************\n")
fout.write( "CPU solve t actual solve t f0 ff ff/f0 Vin\n")
fout.write( " %.2f %.2f %e %e %.5f %.2f\n"%(opt1.solve_t, opt1.wsolve_t,f0, ff, ff/f0, Vin))
fout.write("***************************************************************\n")
fout.write("x0 xf\n")
for i in range(len(x0)):
fout.write("%.15f %.15f\n"%(x0[i],xf[i]))
import shutil
shutil.copy(fi,where)
shutil.copy(fc,where)
In [21]:
print "CPU solve t actual solve t ff/f0"
print " %.2f %.2f %.5f"%(opt1.solve_t, opt1.wsolve_t,ff/f0)
In [22]:
import pickle
pickle.dump(Q00, open(where+"Q0.p", "wb"))
pickle.dump(Q01, open(where+"Qf.p", "wb"))
In [23]:
print where
In [24]:
Vin2 = np.trapz(Q01,x=t)
print Vin2
print Vin
print np.trapz(Q00,x=t)
In [25]:
print "duh"