In [ ]:
%pylab ipympl
Use inline backend followed by jupyter nbconvert rprof_demo.ipynb on command line to create static html to document snap shot.
In [ ]:
# %pylab inline
In [ ]:
import sys, os
sys.path.insert(1, '/user/scratch14_csa/fherwig/repos/PyPPM')
from ppmpy import ppm
from matplotlib import pyplot as plt
import numpy as np
import logging
import math
logging.getLogger("matplotlib").setLevel(logging.ERROR)
In [ ]:
def list_columns(obj, cols=6, columnwise=True, gap=4):
"""
Print the given list in evenly-spaced columns.
Parameters
----------
obj : list
The list to be printed.
cols : int
The number of columns in which the list should be printed.
columnwise : bool, default=True
If True, the items in the list will be printed column-wise.
If False the items in the list will be printed row-wise.
gap : int
The number of spaces that should separate the longest column
item/s from the next column. This is the effective spacing
between columns based on the maximum len() of the list items.
"""
sobj = [str(item) for item in obj]
if cols > len(sobj): cols = len(sobj)
max_len = max([len(item) for item in sobj])
if columnwise: cols = int(math.ceil(float(len(sobj)) / float(cols)))
plist = [sobj[i: i+cols] for i in range(0, len(sobj), cols)]
if columnwise:
if not len(plist[-1]) == cols:
plist[-1].extend(['']*(len(sobj) - len(plist[-1])))
plist = zip(*plist)
printer = '\n'.join([
''.join([c.ljust(max_len + gap) for c in p])
for p in plist])
print(printer)
In [ ]:
dir_repo = '/data/ASDR/PPMstar'
dir_project = 'H-core-M25'
run = 'M35-1536'
path = os.path.join(dir_repo,dir_project,run,'prfs')
rp_set = ppm.RprofSet(path)
In [ ]:
rp_hst = rp_set.get_history()
In [ ]:
rp_hst.get_variables()
In [ ]:
list_columns(rp_hst.get('NDump')[:56], cols=12)
In [ ]:
ifig=1; plt.close(ifig); plt.figure(ifig)
rp_hst.plot_wct_per_dump()
plt.ylim(300, 400)
In [ ]:
rp = rp_set.get_dump(375)
In [ ]:
list_columns(rp.get_header_variables(), cols=3)
In [ ]:
print('Nx = ' + str(rp.get('Nx')))
print('airmu = ' + str(rp.get('airmu')))
print('cldmu = ' + str(rp.get('cldmu')))
Getting header values from the RProfSet instance:
In [ ]:
rp_set.get('Nx',0)
In [ ]:
dump = 375
t=rp_set.get('t',fname=dump)
list_columns(rp_set.get('A', t, num_type='t', resolution='l')[:34], cols=6)
In [ ]:
list_columns(rp_set.get('A', dump, resolution='l')[:34], cols=6)
The argument resolution='l' indicates that we want the low or normal resolution quantity. FV and some others are available in double resolution ('h') which has to do with the high-accuracy advection scheme (cf. Woodward+ 15, Appendix).
These are the rprof column variables available in high and low resolution:
In [ ]:
rp.get_hr_variables()
In [ ]:
list_columns(rp.get_lr_variables(),cols=3)
In [ ]:
# for nice colors, linestyles, etc.
import nugridpy.utils as utils
cb = utils.linestylecb
In the first example we are getting the data and then constructing the plot by hand.
In the second example we use the rp_plot method to make the same plot. It is a bit more concise and convenient, and can ensure that plots have the same look throughout the project.
We are also going to use the rprof method bound_rad to find the radius of the convective boundaruy. Of course, all methods such as get, bound_rad etc have doc strings that can be viewed like this: rp_set.bound_rad?
In [ ]:
ifig=2; plt.close(ifig); plt.figure(ifig)
# get quantities, use high resolution!
fv = rp_set.get('FV',dump,resolution='h')
r = rp_set.get('R',dump,resolution='h')
# plot FV vs. R, semilog it
plt.semilogy(r, fv, color=cb(0)[2], linestyle=cb(0)[0], label='Dump '+str(dump))
# other details in a plot
plt.xlabel('R')
plt.ylabel('FV')
plt.legend()
In [ ]:
rp_set.rp_plot?
In [ ]:
ifig =3 ; plt.close(ifig); plt.figure(ifig)
rp_set.rp_plot([200,400],'FV',ifig=ifig,logy=True)
If using %pylab ipympl the following will plot into the above figure.
In [ ]:
plt.figure(3)
# where is the convective boundary?
conv_bound = rp_set.bound_rad(200,1000,2000,var='FV',criterion='max_grad')
ymin = -3
ymax = 0
plt.vlines(conv_bound[0],ymin,ymax,color='k',\
lw = 0.75,linestyle='--',label='Convective Boundary')
plt.legend(loc=0)
In [ ]:
rp_set.get_computable_quantities()
In [ ]:
N2 = rp_set.compute('N2',dump)
Then you can continue to work with this quantity. The rp_plot method will take computable quantities to plot directly:
In [ ]:
ifig =4 ; plt.close(ifig); plt.figure(ifig)
rp_set.rp_plot(dump,'N2',ifig=ifig)
plt.tight_layout()
In [ ]:
ifig =5 ; plt.close(ifig); plt.figure(ifig)
rp_set.rp_plot(dump,'A',ifig=ifig)
plt.tight_layout()
In [ ]:
1.e6*sqrt(6.e-8)
In [ ]:
rp_set.rprofgui(ifig=19)
In [ ]:
# rp_set.bound_rad?
In [ ]:
Hp=rp_set.compute('Hp',dump)
R = rp_set.get('R',dump)
Ur = rp_set.compute('|Ur|',dump)
rtop = rp_set.bound_rad(dump,1400.,1600.)
In [ ]:
# find convection zone
inds_conv = np.where( R < rtop)[0]
mlen = Hp[inds_conv]
inds=np.where(rtop-R[inds_conv] < Hp[inds_conv])[0]
mlen[inds] = rtop-R[inds_conv][inds]
In [ ]:
alpha_MLT = 1.6
D = (1./3)* alpha_MLT * mlen * Ur[inds_conv] *1.e16 # Mm**2/s
In [ ]:
plt.close(119); plt.figure(119)
plt.plot(R[inds_conv],np.log10(D))
plt.xlabel('$R/\mathrm{Mm}$');plt.ylabel('$\log D / \mathrm{[cm^2/s]}$')
In [ ]:
rp_set.plot_vrad_prof([dump-200,dump],ifig=102)
In [ ]:
dump_min,dump_max,sparse=None,None,10
var = '|U|'
rmin,rmax = 100,150
ylim1,ylim2 = None,None
logy = False
save_pdf = False
tvar = 'time(mins)'
NDump =rp_hst.get('NDump')
dumps = NDump[dump_min:dump_max:sparse]
In [ ]:
timemins = rp_hst.get(tvar)
print("Star time in minutes between dumps: {:.2f}".format(np.average(np.diff(timemins)\
[np.where(np.diff(timemins)>0)])))
# take average of variable between two radii and plot as function of time
things = []
times = []
R = rp_set.get('R', NDump[0])
for dump in dumps:
thing = rp_set.get(var, dump)
times.append(rp_set.get('t',fname=dump)/60.)
things.append(np.average(thing[(R>rmin) * (R<rmax)]))
vars = np.array(things)
times =np.array(times)
i=1
if logy: vars = log10(vars)
plt.close(1093);plt.figure(1093)
plt.plot(times,vars,utils.linestylecb(i)[0],color=utils.colourblind(i),label=None)
plt.ylabel(var);plt.xlabel(tvar)
plt.title("time-"+var+", average "+str(rmin)+" < R/Mm < "+str(rmax))
plt.legend()
plt.close(1092);plt.figure(1092)
plt.plot(dumps,vars,utils.linestylecb(i)[0],color=utils.colourblind(i),label=None)
plt.ylabel(var);plt.xlabel('Ndump')
plt.title("NDump-"+var+", average "+str(rmin)+" < R/Mm < "+str(rmax))
plt.legend()
if save_pdf:
figure(1093)
ylim(ylim1,ylim2)
tight_layout()
savefig("time-"+var+".pdf")
figure(1092)
ylim(ylim1,ylim2)
tight_layout()
savefig("NDump-"+var+".pdf")