In [ ]:
import numpy as np
import os
import fnmatch
import matplotlib
matplotlib.rcParams['figure.dpi'] = 100
matplotlib.rcParams['figure.figsize'] = (12.0, 8.0)
matplotlib.rcParams["font.size"] = 16
import matplotlib.pyplot as plt
from logfile_parse import Logfile
In [ ]:
fnames = ['/d/d9/ychen/2018_production_runs/20180801_L430_rc10_beta07/Group_L430.log',\
'/d/d9/ychen/2018_production_runs/20180802_L438_rc10_beta07/Group_L438.log',\
'/d/d9/ychen/2018_production_runs/20180803_L446_rc10_beta07/Group_L446.log']
logs = {}
for fname in fnames:
logs[fname] = Logfile(fname)
#results = dview.apply_async(Logfile, fname)
for log in logs.values():
log.parse_logfile()
#results = pool.map_async(parse_logfile, log.values())
In [ ]:
for fname, log in logs.items():
fig = plt.figure(figsize=(15,8))
ax = plt.subplot(111)
ax.plot(log.step, log.step_time, '.', markersize=2, c='orange', label='computing time (s) per step')
ax.plot(np.array(log.step), log.tot_blks/10000, '.', markersize=1, label='total # of blocks/10,000')
ax.plot(np.array(log.step), log.tot_leaf_blks/10000, '.', markersize=1, label='total # of leaf blocks/10,000')
ax.vlines(log.startmarks, *ax.get_ylim(), linestyles=':', alpha=0.5)
#for startmark, SU in zip(log.startmarks, log.SUs):
# ax.text(startmark, ax.get_ylim()[1]*0.9, '%.1f SU' % SU)
ax.text(0.05, 0.7, '%.2f SU' % sum(log.SUs), transform=ax.transAxes)
ax.text(0.62, 0.05, log.last_updated.strftime('Last Updated: %Y-%m-%d %H:%M:%S'),
transform=ax.transAxes)
ax.set_ylim(0, 15)
#ax2.set_ylabel('time (sec) /step')
ax.set_xlim(0, log.step[-1])
#ax.set_xlabel('# of step')
#plt.semilogy()
ax.legend(loc=2, numpoints=5, markerscale=4)
#ax2.legend(loc=1, numpoints=5, markerscale=4)
ax.set_title(fname.split('/')[-2])
In [ ]:
for log in logs.values():
repeated_steps = 0
for startmark, endmark in zip(log.startmarks[1:], log.endmarks[:-1]):
repeated_steps += endmark - startmark + 1
print(repeated_steps)
In [ ]:
for fname, log in logs.items():
fig = plt.figure(figsize=(15,8))
ax = plt.subplot(111)
perf = log.tot_blks/log.step_time*3600/log.nProc*48/1E6
ax.plot(log.step, perf, '.', markersize=2, c='green', label='performance')
ax.annotate('%.2f' % np.mean(perf),
(log.step[-len(log.step)//8], np.mean(perf[-len(log.step)//8:])+2))
ax.set_ylim(0, 9)
ax.set_ylabel('# Mblk-step/SU')
ax.set_title(fname.split('/')[-2])
ax2 = ax.twinx()
ax2.fill_between(log.step, log.min_blks, log.max_blks, label='# of blocks / proc')
ax2.fill_between(log.step, log.min_leaf_blks, log.max_leaf_blks, label='# of leaf blocks / proc')
ax2.set_xlim(0, log.step[-1])
ax2.set_ylim(0, 600)
ax2.vlines(log.startmarks, *ax2.get_ylim(), linestyles=':', alpha=0.5)
ax.legend(loc=2, numpoints=5, markerscale=2)
ax2.legend(loc=1)
In [ ]:
from scipy import optimize
def powerlaw(xx, a, p):
return a*xx**p
def fit_powerlaw(times, cmulativeSUs):
popt, pcov = optimize.curve_fit(powerlaw, times, cumulativeSUs)
return popt
def SUs_proj(fname, popt):
if 'L430' in fname:
times_proj = np.linspace(0, 25, 100)
elif 'L438' in fname:
times_proj = np.linspace(0, 65, 100)
elif 'L446' in fname:
times_proj = np.linspace(0, 75, 100)
print(fname, popt)
#a, p = popt
SUs = powerlaw(times_proj, *popt)
return times_proj, SUs
labels = {'20180801_L430_rc10_beta07': 'low power',
'20180802_L438_rc10_beta07': 'mid power',
'20180803_L446_rc10_beta07': 'high power',}
for fname, log in logs.items():
times = [0]
cumulativeSUs = [0]
SUs = 0
for start, end, SU in zip(log.startmarks, log.endmarks, log.SUs):
i = np.argwhere(log.step == start)[0,0]
j = np.argwhere(log.step == end)[0,0]
#print('n_start = %6i(%5.2f Myr) , n_end = %6i(%5.2f Myr) : %5.2f SUs' %
# (start, log.t_sim[i]/3.15576E13, end, log.t_sim[j]/3.15576E13, SU))
times.append(log.t[j]/3.15576E13)
SUs += SU
cumulativeSUs.append(SUs)
popt = fit_powerlaw(times, cumulativeSUs)
print(popt)
#SUs_proj = 5*times_proj**1.8
#SUs_proj = 1.1*times_proj**1.8
times_proj, SUs = SUs_proj(fname, popt)
label = labels[fname.split('/')[-2]] + r' - %5i SUs'
line, = plt.plot(times, cumulativeSUs, 'o-', label=label % cumulativeSUs[-1])
plt.plot(times_proj, SUs, ':', c=line.get_color(),\
label=r'(Projected) %.2f$\times t^{%.2f}$ - %6i SUs' % (tuple(popt)+(SUs[-1],)))
plt.xlabel('sim time (Myr)')
plt.ylabel('SUs')
plt.legend(fontsize=12)
#plt.xlim(0,20)
#plt.ylim(0,1000)
In [ ]:
def get_nozzledata(dirname):
try:
nozzledata = np.genfromtxt(os.path.join(dirname, 'nozzleVec.dat'))
ttnoz = nozzledata[:,1]
xx = nozzledata[:,2]
yy = nozzledata[:,3]
zz = nozzledata[:,4]
thetas = np.arccos(zz)/np.pi*180
phis = np.arctan2(yy, xx)
return ttnoz, xx, yy, zz, thetas, phis
except:
pass
def get_dat(dirname):
fdat = fnmatch.filter(os.listdir(dirname), 'Group_L???.dat')[0]
data = np.genfromtxt(os.path.join(dirname, fdat))
tt = data[:,0]
Masses = data[:,1]
Etots = data[:,5]
Ekins = data[:,6]
Eints = data[:,7]
Emags = data[:,8]
Ptots = (Etots[1:]-Etots[:-1])/(tt[1:]-tt[:-1])
Pkins = (Ekins[1:]-Ekins[:-1])/(tt[1:]-tt[:-1])
Pints = (Eints[1:]-Eints[:-1])/(tt[1:]-tt[:-1])
#print max(Ptots), max(Pkins), max(Pints)
return tt, Masses, Etots, Ekins, Eints, Emags, Ptots, Pkins, Pints
In [ ]:
dirnames = ['/d/d9/ychen/2018_production_runs/20180801_L430_rc10_beta07/',\
'/d/d9/ychen/2018_production_runs/20180802_L438_rc10_beta07/',\
'/d/d9/ychen/2018_production_runs/20180803_L446_rc10_beta07/']
table = []
for dirname in dirnames:
table.append(get_dat(dirname))
In [ ]:
names = ['Low power', 'Mid power', 'High power']
#colors = ['yellow', 'r', 'b', 'g', 'cyan']
plt.figure(figsize=(10,8))
for i in range(len(table)):
tt, Masses, Etots, Ekins, Eints, Emags, Ptots, Pkins, Pints = table[i]
#plt.plot(tt[1:]/3.154E13, (tt[1:]-tt[:-1]), label=names[i], lw=0, marker='.', markersize=1, alpha=0.7)
plt.plot(tt[1:]/3.154E13, (tt[1:]-tt[:-1]), label=names[i], lw=0, marker='.', markersize=2, alpha=0.7)
#plt.scatter(tt[:], Etots[:], s=1, linewidth=0, label='dt')
#plt.ylim(2.563E62, 2.563E62+9E59)
#plt.ylim(0, 4E58)
plt.ylim(4E7,2E10)
#plt.xlim(0, 2)
plt.semilogy()
plt.ylabel('dt (s)')
plt.xlabel('t (Myr)')
plt.legend(loc=4, fontsize=12, numpoints=5, markerscale=2)
ymin, ymax = plt.gca().get_ylim()
ax2 = plt.twinx()
ax2.set_ylim(ymin/3.15E7, ymax/3.15E7)
ax2.semilogy()
ax2.set_ylabel('dt (yr)')
plt.show()
In [ ]:
i=2
L = 4E44
tmax = 40
tt, Masses, Etots, Ekins, Eints, Emags, Ptots, Pkins, Pints = table[i]
dirname = dirnames[i]
fig = plt.figure(figsize=(8,6))
plt.scatter(tt[1:]/3.15569E13, (Etots[1:]-Etots[:-1])/(tt[1:]-tt[:-1]), s=1, c='b', linewidth=0, label='Total Power')
plt.scatter(tt[1:]/3.15569E13, (Ekins[1:]-Ekins[:-1])/(tt[1:]-tt[:-1]), s=1, c='r', linewidth=0, label='Kinetic Power')
plt.scatter(tt[1:]/3.15569E13, (Eints[1:]-Eints[:-1])/(tt[1:]-tt[:-1]), s=1, c='g', linewidth=0, label='Internal Power')
plt.scatter(tt[1:]/3.15569E13, (Emags[1:]-Emags[:-1])/(tt[1:]-tt[:-1]), s=1, c='purple', linewidth=0, label='Magnetic Power')
#plt.ylim(1E44, 1.5E45)
plt.ylim(-1E43, L*1.25)
plt.ylabel('power (erg/s)')
plt.axhline(0, c='grey', alpha=0.5)
plt.axhline(L, lw=0.5, c='grey')
plt.axhline(L*1.1, lw=0.5, c='grey')
plt.legend(loc=2)
#plt.semilogy()
#plt.twinx()
#plt.scatter(ttnoz/3.15569E13, thetas, s=1, linewidth=0, c='y', alpha=0.3, label='theta')
#plt.scatter(ttnoz, phis, s=1, linewidth=0, c='g', alpha=0.5, label='phi')
plt.legend(loc=1, numpoints=5, markerscale=4)
plt.xlim(0, tmax)
plt.xlabel('t (Myr)')
plt.title(dirname.split('/')[-1])
plt.show()
In [ ]: