Data analysis demo for H-core-M25 stellar hydro project

Last update: Nov 25, 2019.

This notebook contains a demonstration how to analyse the 3D filtered moms data with the 1D radial profile rprof data from PPMstar 3D hydrodynamic simulations.

Data for this demo

The examples are for the project H-core-M25 (this is the project identifier), the H-core convection simulations of a 25M$_\odot$ star.

Two runs are used:

M29: $768^3$ grid
M35: $1536^3$ grid

M29, M35 are the run identifier. Keep run and project identifier attached to all derived data products.

Both runs have 1000x heating which increases their convective velocities by a factor of 10.

For each run there are two types of data to be read for this demo:

moms data is the spatially filtered data (2-byte data on reduced grid by factor four in each direction) in 3D
rprofs data are spherically averaged radial profiles

Location of data

The data is staged on the UVic Astrophysics Simulation Data Repository (ASDR) mounted in /data/ASDR. The repository contains the project folder H-core-M25.

Python asumptions

The server defaults each notebook to %pylab ipympl



In [ ]:

    
## use this for final run to export with images to pdf, markdown or html
#%pylab inline



In [ ]:

    
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as color
import nugridpy.utils as utils
import logging
import collections
import itertools
import sys
import os

# you can use sys.path.insert to add your own cloned copy of PyPPM which you can modify
# sys.path.insert(0,'/home/user/home/PyPPM/')
from ppmpy import ppm

# plotting utilities
ifig = 0
ptrack = {}
cb = utils.linestylecb # colours

# create smaller labelsizes on axes
plt.rc('xtick', labelsize=7)
plt.rc('ytick', labelsize=7)

# named tuple for using rprofs and momsdata
hydro = collections.namedtuple('hydro', ['moms','rprof'])

# turn off matplotlib messages
logging.getLogger("matplotlib").setLevel(logging.CRITICAL)



In [ ]:

    
%%bash
ls /data/ASDR/PPMstar/H-core-M25/

What quantities have what index?

Due to how the moments data are decompressed the quantities written into the moms data file can be written in any order and do not contain any metadata to notify the MomsDataSet class what quantities are where in the large array. The quantities can be called with an index or their name IF you define the $\textit{var_list}$ variable in the MomsDataSet instantiation. For the runs being looked at here, the momsdata indices and their appropriate quantities are:

index	quantity
0	x
1	$\vec{u_{x}}$
2	$\vec{u_{y}}$
3	$\vec{u_{z}}$
4	$\lvert\vec{u_{t}}\rvert$
5	$\lvert\vec{u_{r}}\rvert$
6	$\lvert\vec{\omega}\rvert$
7	P
8	rho
9	fv

Note that these are just 10 out of 32 quantities that can be made available in the moms data.
fv is the fractional volume of the material initially only outside the convection zone.

Some Helpful Definitions

$\mu$ = fv $\times$ $\mu_{\mathrm{cld}}$ + (1 - fv) $\times$ $\mu_{\mathrm{air}}$

T = $\frac{P \mu}{\rho R_{gas}}$

$R_{gas}$ = 8.314462

$\vec{\omega}$ = $\vec{\nabla} \times \vec{u}$

$R_{gas}$, and all other quantities, are in code units which are defined in this file:



In [ ]:

    
%%bash
cat ./code-units.txt

We can find out how to instantiate the MomsDataSet object within the notebook as:



In [ ]:

    
ppm.MomsDataSet?

All class methods contain doc strings that describe the arguments and outputs.



In [ ]:

    
# paths to data
dir_repo    = '/data/ASDR/PPMstar'
dir_project = 'H-core-M25'

# the runs and their resolutions
runs = ['M29', 'M35']
resolution = ['768', '1536']

# dumps for each run we could analyze
dumps = [650, 375]

# are we using high resolution?
use_highres = False

if use_highres:
    use_run = [False, True]
    
else:
    use_run = [True, False]

# grab the appropriate runid, runresolution and dump
runid, runresolution, mydump = itertools.compress(runs+resolution+dumps, use_run*3)

# create the hydro object
var_list = ['xc','ux','uy','uz','|ur|','|ut|','|w|','P','rho','fv']
dir_repo = os.path.join(dir_repo, dir_project, '{:s}-{:s}'.format(runid, runresolution))
                        
myhydro = hydro(ppm.MomsDataSet(dir_name=os.path.join(dir_repo, 'myavsbq'), init_dump_read=mydump, 
                                dumps_in_mem=2, var_list=var_list, 
                                rprofset=ppm.RprofSet(os.path.join(dir_repo, 'prfs'))),
                ppm.RprofSet(os.path.join(dir_repo, 'prfs')))

Basic grid properties



In [ ]:

    
# 3D datacube of the cartesian coordinates at this dump
x,y,z=myhydro.moms.get_cgrid()

# it is formatted such that var[z,y,x] refers to the variable at the (z,y,x) coordinates that those 
# indices correspond to. The indexing goes from negative to positive 
moms_ngridpoints = myhydro.moms.moms_ngridpoints
print('The unique x,y,z coordinates are:')
print(x[np.random.randint(0,moms_ngridpoints), np.random.randint(0,moms_ngridpoints), :])

# to convert to spherical coordinates
r,theta,phi=myhydro.moms.get_sgrid()

Grid and time resolution of moms data and rprofs (PPMStar simulation)



In [ ]:

    
# breakdown of cells and shapes of the simulations
ppmstar_dim = len(myhydro.rprof.get('R',fname=0,resolution='h'))
print('The original PPMStar simulation was run with a cubic dimension of ('+','.join(3*[str(ppmstar_dim)])+
     ') while the moments data is downsampled by a factor of 4 with a cubic dimension of ('+
      ','.join(map(str,x.shape))+')')

print('')

# spatial resolution for momsdata compared with rprofs (simulation)
deex_ppmstar = myhydro.rprof.get('deex',fname=0)
print('The grid resolution of the original PPMStar simulation is {:.2f}Mm while the moms data is {:.2f}Mm'
      .format(deex_ppmstar,np.diff(x[0,0,:])[0]))



In [ ]:

    
# temporal resolution is in the rprof data
print('The temporal resolution of the moms data is the same as the PPMStar output which averages around {:0.2f}'.
      format(np.mean(np.diff(myhydro.rprof.get_history().get('time(mins)')))),'minutes of star time per dump')

print('')

print('The run-time temporal resolution of the PPMStar output averages around {:0.2f}'.
      format(np.mean(myhydro.rprof.get_history().get('dt(secs)'))),'seconds per cycle')

Histogram of radii

increasing to 1/2 length of grid, then decreasing as only fraction of shell in box



In [ ]:

    
# this ensures that I can open an close the figure that was plotted within this cell. 
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 0

# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
    plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
else:
    ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
    
# set the colour/linestyle counter to 0 and create axes
cbc = 0
ax = fig.add_subplot(111)

# histogram of radius values
hist_data = ax.hist(r.flatten(),86)

# plot details
xmax = np.max(x)
ylims = [0,np.max(hist_data[0])]
ax.vlines(xmax,*ylims,linestyles='--',lw=0.5)
ax.set_xlabel('r / Mm')
ax.set_ylabel('')
ax.set_ylim(ylims)

Planar slice image



In [ ]:

    
# extent in x,y
extent=[np.min(x),np.max(x),np.min(y),np.max(y)]

# the take a slice at z=0. This is halfway through the indices
slice_num = int(moms_ngridpoints/2)

Radius



In [ ]:

    
# this ensures that I can open an close the figure that was plotted within this cell. 
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 1

# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
    plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
else:
    ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
    
# set the colour/linestyle counter to 0 and create axes
cbc = 0
ax = plt.axes([0.01, 0.98, 0.88, 0.98])

ax.imshow(r[:,:,slice_num],extent=extent)

# plot details
ax.set_ylabel('y / Mm')
ax.set_xlabel('x / Mm')

# create a colorbar for the imshow data, use a new axes for it
cax = fig.add_axes([0.90, 0.98, 0.08, 0.98])
cbar1 = fig.colorbar(ax.images[0], cax=cax, orientation='vertical',format='%d',
                     label=r'r / Mm')

$| u_{\perp} |$

We can simply construct the magnitude of the tangential velocity from the cartesian vector velocity components



In [ ]:

    
# get the cartesian velocity components
ux = myhydro.moms.get('ux',fname=mydump)
uy = myhydro.moms.get('uy',fname=mydump)
uz = myhydro.moms.get('uz',fname=mydump)

# convert to spherical components
ur, utheta, uphi = myhydro.moms.get_spherical_components(ux, uy, uz)
u_perp = np.sqrt(np.power(utheta,2.0) + np.power(uphi,2.0))

# store for later, |u|
u_cartesian = myhydro.moms.norm(ux, uy, uz)
u_spherical = myhydro.moms.norm(ur, utheta, uphi)

# they are practically the same
print((u_cartesian - u_spherical).mean(), (u_cartesian - u_spherical).std())



In [ ]:

    
# this ensures that I can open an close the figure that was plotted within this cell. 
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 2

# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
    plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
else:
    ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
    
# set the colour/linestyle counter to 0 and create axes
cbc = 0
ax = plt.axes([0.01, 0.98, 0.88, 0.98])

# get the tangential velocity and make it km/s rather than Mm/s
ax.imshow(u_perp[:,:,slice_num]*1e3,extent=extent,vmin=0)

# plot details
ax.set_ylabel('y / Mm')
ax.set_xlabel('x / Mm')

# create a colorbar for the imshow data, use a new axes for it
cax = fig.add_axes([0.90, 0.98, 0.08, 0.98])
cbar1 = fig.colorbar(ax.images[0], cax=cax, orientation='vertical',format='%.0f',
                     label=r'$u_{\perp}$ / km s$^{-1}$')

Radial profiles

Radial profiles can be directly taken from the rprof data sets. They can also be constructed from the moms data. This is demonstrated below.



In [ ]:

    
# for rprofs we will grab the spherical averages directly
T9_rprof = myhydro.rprof.compute_T9(fname=mydump)
R_rprof = myhydro.rprof.get('R',fname=mydump,resolution='l')

# for moments data we will have to use the base components and compute T9
cldmu = myhydro.rprof.get('cldmu',fname=0)
airmu = myhydro.rprof.get('airmu',fname=0)
Rgas = 8.314462

# calculate mu for T
fv = myhydro.moms.get('fv',fname=mydump)
mu = cldmu * fv + (1-fv)*airmu

P = myhydro.moms.get('P',fname=mydump)
rho = myhydro.moms.get('rho',fname=mydump)
T9 = (mu * P) / (Rgas * rho)

# get the radial profile of T
T9_moms, R_moms = myhydro.moms.get_rprof(T9)

T9



In [ ]:

    
# this ensures that I can open an close the figure that was plotted within this cell. 
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 3

# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
    plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
else:
    ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
    
# set the colour/linestyle counter to 0 and create axes
cbc = 0
ax = fig.add_subplot(111)

ax.plot(R_rprof, T9_rprof, label='Rprof', color=cb(cbc)[2], ls=cb(cbc)[0])
cbc += 1

ax.plot(R_moms, T9_moms, label='MomsData', color=cb(cbc)[2], ls=cb(cbc)[0])
cbc += 1

# plot details
ax.set_ylabel('T9')
ax.set_xlabel('r / Mm')
ax.set_xlim([0, R_rprof.max()])

ax.legend(frameon=False,fontsize=8)

|u|



In [ ]:

    
U_rprof = myhydro.rprof.get('|U|',fname=mydump)

# we computed |U| from the cartesian components of U as well as spherical and they are the same on the 
# entire grid. Use cartesian
U_moms, R_moms = myhydro.moms.get_rprof(u_cartesian)



In [ ]:

    
# this ensures that I can open an close the figure that was plotted within this cell. 
# It doesn't matter in what order I run the notebook cells, this keeps track of it
celln = 4

# configure plots
close_local = ppm.close_plot(celln,ifig,ptrack)
if close_local[0]:
    plt.close(fig); ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
else:
    ifig += 1; fig = plt.figure(ifig,figsize=(3.5,2.5),dpi=300)
    ppm.add_plot(celln,ifig,ptrack)
    
# set the colour/linestyle counter to 0 and create axes
cbc = 0
ax = fig.add_subplot(111)

# convert to km/s
ax.plot(R_rprof, U_rprof*1e3, label='Rprof', color=cb(cbc)[2], ls=cb(cbc)[0])
cbc += 1

ax.plot(R_moms, U_moms*1e3, label='MomsData', color=cb(cbc)[2], ls=cb(cbc)[0])
cbc += 1

# plot details
ax.set_ylabel('|U| / km s$^{-1}$')
ax.set_xlabel('r / Mm')
ax.set_xlim([0, R_rprof.max()])
ax.set_ylim([0,16])

ax.legend(frameon=False,fontsize=8)



In [ ]: