In [11]:
%load_ext autoreload
%autoreload 2
import sys
sys.path.append('/user/andrassy/PyPPM')
from ppmpy import ppm
from matplotlib import pyplot as plt
%matplotlib nbagg
from nugridpy import utils
cb = utils.colourblind
In [12]:
yprof_paths = {'D2':'/data/ppm_rpod2/YProfiles/O-shell-M25/D2/', \
'E10':'/data/ppm_rpod2/YProfiles/RAWD/E10/', \
'F4':'/data/ppm_rpod2/YProfiles/AGBTP_M2.0Z1.e-5/F4/'}
yprofs = {rid:ppm.yprofile(path) for rid, path in yprof_paths.items()}
rprof_paths = {'M12':'/user/niagara_projects/PPM2.0/M_H-core-M25/M12-Hcore25-1536-1000x/prfs/', \
'N3':'/user/niagara_projects/PPM2.0/N_LowZRARD/N3-LowZRAWD-1152-10x/prfs/'}
rp_sets = {rid:ppm.RprofSet(path) for rid, path in rprof_paths.items()}
run_data = {**yprofs, **rp_sets}
This method can find the boundary radius as defined by a local minimum or maximum in the gradient of a variable or by searching for a certain value of a variable (e.g. FV = 0.5). The method can also return the scale height of that variable at the boundary, which entrainment_rate() makes use of.
In [13]:
rid = 'D2'
dumps = range(84, 164, 1)
r_min = 7.5
r_max = 8.5
var = 'ut'
criterion = 'min_grad'
rb, H = run_data[rid].boundary_radius(dumps, r_min, r_max, var=var, \
criterion=criterion,
return_var_scale_height=True)
offset = -1.
rt = rb + offset*np.abs(H)
ifig=1; plt.close(ifig); plt.figure(ifig)
plt.plot(dumps, rt, '-', color=cb(5), label=r'r$_\mathrm{top}$')
plt.plot(dumps, rb, '--', color=cb(8), label=r'r$_\mathrm{b}$')
plt.legend(loc=0)
plt.xlabel('dump #')
plt.ylabel('r [Mm]')
plt.title(rid)
Out[13]:
In [14]:
rid = 'M12'
dumps = range(200, 492, 5)
r_min = 1450.
r_max = 1600.
var = 'FV'
criterion = 'max_grad'
rb, H = run_data[rid].boundary_radius(dumps, r_min, r_max, var=var, \
criterion=criterion,
return_var_scale_height=True)
offset = -1.
rt = rb + offset*np.abs(H)
ifig=2; plt.close(ifig); plt.figure(ifig)
plt.plot(dumps, rt, '-', color=cb(5), label=r'r$_\mathrm{top}$')
plt.plot(dumps, rb, '--', color=cb(8), label=r'r$_\mathrm{b}$')
plt.legend(loc=0)
plt.xlabel('dump #')
plt.ylabel('r [Mm]')
plt.title(rid)
Out[14]:
In [15]:
rid = 'D2'
dumps = range(84, 164, 1)
r_min = 7.5
r_max = 8.5
var = 'ut'
criterion = 'min_grad'
offset = -1.
_ = run_data[rid].entrainment_rate(dumps, r_min, r_max, var=var, \
criterion=criterion, \
offset=offset, ifig0=3)
In [16]:
rid = 'M12'
dumps = range(200, 492, 5)
r_min = 1450.
r_max = 1600.
var = 'FV'
criterion = 'max_grad'
offset = -1.
_ = run_data[rid].entrainment_rate(dumps, r_min, r_max, var=var, \
criterion=criterion, \
offset=offset, ifig0=5)
Run N3 includes H burning with a temperature correction. The user has to specify burn_func, which is a function that computes the mass burning rate per unit volume as a function of radius. The first argument of this function must be the cycle (dump) number and the remaining parameters are specified in burn_args.
In [17]:
rid = 'N3'
dumps = range(0, 400, 5)
r_min = 22.
r_max = 26.
var = 'FV'
criterion = 'max_grad'
offset = -2.
fit_interval = 60.*np.array([200., 370.])
burn_func = run_data[rid].compute_rhodot_C12pg
T9corr_params = {'kind':1, 'params':{'a':0.46, 'b':0.77}}
burn_args = {'T9corr_params':T9corr_params}
res = run_data[rid].entrainment_rate(dumps, r_min, r_max, var=var, \
criterion=criterion, \
burn_func=burn_func, \
burn_args=burn_args, offset=offset, \
fit_interval=fit_interval, ifig0=7)
In [18]:
rid = 'E10'
dumps = range(150, 427, 5)
r_min = 32.
r_max = 42.
var = 'ut'
criterion = 'min_grad'
offset = -1.
fit_interval = 60.*np.array([220., 400.])
burn_func = run_data[rid].compute_rhodot_C12pg
T9corr_params = {'kind':1, 'params':{'a':0.6, 'b':0.85}}
burn_args = {'T9corr_params':T9corr_params, \
'airmu':1.4, \
'cldmu':0.3, \
'fkair':0.203606102635, \
'fkcld':0.885906040268, \
'atomicnoair':6.65742024965, \
'atomicnocld':1.34228187919}
res = run_data[rid].entrainment_rate(dumps, r_min, r_max, var=var, \
criterion=criterion, \
burn_func=burn_func, \
burn_args=burn_args, offset=offset, \
fit_interval=fit_interval, ifig0=9)
F4 is a more complicated case, because H burning in this run was off in the first 616 dumps and on afterwards. We define a custom burn_func(), which receives the data object as well as all burning parameters via burn_args. In principle, one could define burn_func() that computes the total of two different reactions in the same way.
In [19]:
rid = 'F4'
dumps = range(50, 1350, 10)
r_min = 27.7
r_max = 28.5
var = 'ut'
criterion = 'min_grad'
offset = -1.
fit_interval = 60.*np.array([1100., 1550.])
burn_args = {'data_object':run_data[rid], \
'airmu':1.39165, \
'cldmu':0.725, \
'fkair':0.203606102635, \
'fkcld':0.885906040268, \
'atomicnoair':6.65742024965, \
'atomicnocld':1.34228187919}
def burn_func(cycle, **kwargs):
r = kwargs['data_object'].get('Y', 0, resolution='l')
burn_rate = np.zeros(len(r))
data_object = kwargs['data_object']
burn_args = dict(kwargs)
del burn_args['data_object']
if cycle > 620:
burn_rate = data_object.compute('rhodot_C12pg', cycle, \
extra_args=burn_args)
return burn_rate
res = run_data[rid].entrainment_rate(dumps, r_min, r_max, var=var, \
criterion=criterion, burn_func=burn_func, \
burn_args=burn_args, offset=offset, \
fit_interval=fit_interval, ifig0=11)
In [ ]: