{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline\n", "import yt \n", "yt.mylog.setLevel(50) #suppress yt outputs " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To analyze the hdf5 files that the simulation outputs,we will be using yt, a well-supported Python package for visualizing volumetric, AMR data from astrophysical simulations. For more information, refer to the excellent [yt documentation](http://yt-project.org/docs/3.2/) and [series of notebook tutorials](http://yt-project.org/docs/dev/quickstart/index.html) that is much more comprehensive than what is covered here. In this notebook, we will be covering yt basics, as well as using some custom scripts used for analyzing the FLASH data in ``plotSim.py`` that I wrote." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/anaconda2/lib/python2.7/site-packages/matplotlib/__init__.py:1350: UserWarning: This call to matplotlib.use() has no effect\n", "because the backend has already been chosen;\n", "matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n", "or matplotlib.backends is imported for the first time.\n", "\n", " warnings.warn(_use_error_msg)\n" ] } ], "source": [ "from scripts.plotSim import *" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/global/u2/d/dorislee/data_astroSim/data\n" ] } ], "source": [ "cd ~/data_astroSim/data/" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exploratory Data Analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Loading in the dataset for the last timestep in our Rotating Sink Sphere simulation: " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ds = yt.load(\"sphere_hdf5_chk_0494\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By default, FLASH outputs its data in cgs (except in MHD..), so the code units in yt is actually in cgs." ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0 cm 1.0 g 1.0 s\n" ] } ], "source": [ "print ds.length_unit,ds.mass_unit,ds.time_unit " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Often useful information about a simulation timestep output: \n", " " ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Current simulation timestep: 4.94073212769e+12 code_time\n" ] } ], "source": [ "print \"Current simulation timestep:\", ds.current_time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is a list of variables that are stored in the simulation outputs, note that some variables (e.g mag) will not have data (since this is a hydro run):" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['dens', 'divb', 'eint', 'ener', 'gamc', 'game', 'gpol', 'gpot', 'grac', 'magp', 'magx', 'magy', 'magz', 'pden', 'pres', 'sgax', 'sgay', 'sgaz', 'sgxo', 'sgyo', 'sgzo', 'shok', 'temp', 'velx', 'vely', 'velz']\n" ] } ], "source": [ "print map(lambda x:str(x[1]),filter(lambda x :x[0]=='flash',ds.field_list))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is just a demo of the long list of variables that yt can derive from the few variables stored in the simulation as shown above:" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['particle_position_cylindrical_z', 'cylindrical_theta', 'particle_angular_momentum_z', u'game', 'particle_radius', u'divb', 'total_energy', 'particle_position_cylindrical_radius', 'density_gradient_magnitude', 'particle_cylindrical_velocity_theta', 'pressure', 'particle_velocity_z', 'particle_position_spherical_radius', 'magnetic_field_strength', 'io_nn_velocity_x']\n" ] } ], "source": [ "print map(lambda x: x[1],ds.field_info.keys())[:15]" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# SlicePlots" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We many things we can do on the yt dataset ``ds``, such as making SlicePlots, create a volumetric rendering, or making radial profiles by cutting through the simulation box." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we make a SlicePlot of the density through the z direction" ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slc = yt.SlicePlot(ds, \"z\",\"density\")\n", "slc.set_figure_size(5) #sets the size of the plot\n", "slc.set_cmap(\"all\",\"rainbow\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot what the AMR grid looks like. In this simulation, we are using Jean's criterion based AMR refinement. The AMR structure would look quite different if we used density based refinement (which refines based on density gradients across cells). Note that the grid lines that yt is drawing denotes where the *blocks* are __NOT__ where the *cells* are. (Recall that each block consists of 8 cells)" ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 67, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slc.annotate_grids()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also do fancy things like zooming in and plotting the velocity vectors." ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slc.zoom(4)\n", "slc.annotate_grids()\n", "slc.annotate_velocity(normalize=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "zlim defines the limits on the colorbar for the variable that you are plotting. This is sometimes useful when we are looking at data across many timesteps." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "~~~python\n", "slc.set_zlim(\"density\",zmin,zmax)\n", "\n", "~~~" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_dens(0,fname=\"sphere\")" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_var(0,\"temperature\", \"sphere\")\n", "plot_var(0,\"pressure\", \"sphere\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Without fixing the colorbar," ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_dens(300,fname=\"sphere\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fixing the colorbar," ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_dens(300,fname=\"sphere\",zmin = 5e-21,zmax = 1e-17)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As shown above, fixing the colorbar allows us to compare our simulation against the initial conditions. For example we may be interested in seeing which areas of the sphere experience a density increase compared to the initial sphere, rather than by *how much* does the density increase. On the other hand, allowing the colorbar to vary enable us to see more detailed structures that might be saturated if we retained the initial color scheme. It is important to think about whether the visualization settings may affect what we see in the data. \n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_dens(494,fname=\"sphere\",zoom=8)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_dens(494,fname=\"sphere\",zoom=8,zmin = 5e-21,zmax = 1e-17)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Adding user-defined parameters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Even though yt has a comprehensive list of derived variables , sometimes it is useful to define our own parameters. For example, Foster and Chevalier (1993) defines a set of non-dimensional units used for their plots, so in order to compare our simulation result to theirs, we must additionally convert into these nondimensional units. \n", "\\begin{align}\n", "D = \\frac{\\rho}{\\rho_c} = 9.09\\times 10^{18}\\rho\n", "\\\\ m = \\frac{GM\\sqrt{4\\pi G\\rho_c}}{a^3}=8.54\\times10^{-34} M\n", "\\\\v = \\frac{u}{a} = 3.48\\times10^{-5} u \n", "\\\\ \\xi =\\frac{\\sqrt{4\\pi G \\rho_c}}{a} r = 1.0568\\times10^{-17} r \n", "\\\\ \\tau = \\sqrt{4\\pi G \\rho_c} t = 3.036\\times 10^{-13} t \n", "\\end{align}\n", "where $\\rho$, M, u, r,t are in cgs units.\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline\n", "import yt\n", "yt.mylog.setLevel(50)\n", "from yt.units import second, g, cm ,dyne\n", "G = 6.674e-8*cm**3/second**2/g\n", "rho_c = 1.1e-17*g/cm**3\n", "conversion = sqrt(4*pi*G*rho_c)\n", "aT = 28730*cm/second # for 10K gas \n", "ctr= 5e18*cm\n", "def _radius(field,data):\n", " radius = sqrt((data[\"x\"]-ctr)**2+(data[\"y\"]-ctr)**2+(data[\"z\"]-ctr)**2)\n", "def _norm_r(field,data):\n", " return data[\"radius\"]*conversion/aT\n", "def _norm_v(field,data):\n", " return data[\"radial_velocity\"]/aT\n", "def _norm_d(field,data):\n", " return data[\"density\"]/rho_c\n", "yt.add_field(\"radius\",function= _radius,units=\"cm\")\n", "yt.add_field(\"norm_r\",function= _norm_r,units=\"\")\n", "yt.add_field(\"norm_v\",function= _norm_v,units=\"\")\n", "yt.add_field(\"norm_d\",function= _norm_d,units=\"\")\n", "ds = yt.load(\"../../data_astroSim/data/sphere_hdf5_chk_0000\")\n", "data = ds.all_data()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the units are non-dimensional " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 908.41256519 903.6691589 898.97660951 ..., 898.97660951 903.6691589\n", " 908.41256519] dimensionless\n", "[ 9.89269385e-19 4.94634692e-19 1.48390408e-18 ..., 0.00000000e+00\n", " -4.94634692e-19 0.00000000e+00] dimensionless\n", "[ 0.00045455 0.00045455 0.00045455 ..., 0.00045455 0.00045455\n", " 0.00045455] dimensionless\n" ] } ], "source": [ "print data[\"norm_r\"]\n", "print data[\"norm_v\"]\n", "print data[\"norm_d\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now use these non-dimensional, derived fields in a profile plot: " ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sphere = ds.sphere(\"c\", (0.5, \"pc\"))\n", "prof = yt.ProfilePlot(sphere, \"norm_r\",\"norm_d\")\n", "prof.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# FixedResolutionBuffer" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "turns your AMR data into a 2D or 3D numpy array which you can then use for plotting in \n", "If you are used to matplotlib or other plotting library \n", "I find that generally plotting with --- yt is sufficient and usually better than whatever ----, but FRBs are especially useful for quantitatively analyzing the data. For example in scripts/----, M/R plots " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "Inside, scripts/" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/global/project/projectdirs/astro250/dlee/FLASH4.3/object\n" ] } ], "source": [ "cd ~/proj/dlee/FLASH4.3/object/" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "sphere.dat sphere_hdf5_part_0000 sphere.log\r\n", "sphere_hdf5_chk_0000 sphere_hdf5_plt_cnt_0000\r\n" ] } ], "source": [ "ls sphere*" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_dens(0,\"sphere\")" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_var(0,\"temperature\", \"sphere\")" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_var(0,\"pressure\", \"sphere\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }