Plonk tutorial for the 3rd Phantom + MCFOST users workshop, Monash University, Feb 25, 2020 -- https://phantomsph.bitbucket.io/workshop2020

_README.md

This Plonk tutorial was written for Plonk v0.3.0. Parts may no longer work for later versions of Plonk. Please see the documentation at https://plonk.readthedocs.io/ for current documentation and examples.

installation-instructions.md

Installing Plonk in a new conda environment

In general, I recommend setting conda-forge as the main channel from which to install all packages in all environments. The following commands will do this:

conda config --set channel_priority true
conda config --add channels conda-forge

Now, create a new conda environment to install Plonk in:

conda create --name myenv

The environment name is myenv. Feel free to change it as you like.

Now activate the environment:

conda activate myenv

You can install Plonk in the environment with:

conda install plonk

Don't forget that if you open a new terminal or launch Jupyter lab you must do it after activating the environment in which you installed Plonk.

Plonk-tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plonk tutorial\n",
    "\n",
    "by Daniel Mentiplay ([@dmentipl](https://dmentipl.github.io/))\n",
    "\n",
    "This notebook contains a Plonk tutorial written for the 3rd Phantom + MCFOST users workshop, Monash University, Feb 24-28, 2020. See <https://phantomsph.bitbucket.io/workshop2020> for details on the workshop. The material here is derived from the Plonk documentation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is Plonk?\n",
    "> Plonk is open source Python package for analysis and visualization of smoothed particle hydrodynamics data.\n",
    "\n",
    "- For detailed Plonk usage and installation instructions, see the documentation at <https://plonk.readthedocs.io/>.\n",
    "- The source code is available on GitHub at <https://github.com/dmentipl/plonk>.\n",
    "- The peer-reviewed paper is in JOSS at <https://joss.theoj.org/papers/10.21105/joss.01884>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Supported data formats\n",
    "\n",
    "🚨 **Plonk (v0.3.0, latest as of Feb 2020) can *only* read Phantom snapshots in the HDF5 format.** 🚨\n",
    "\n",
    "There are detailed instructions on running Phantom with HDF5 and converting standard snapshots to HDF5 in the documentation at <https://phantomsph.readthedocs.io/en/latest/hdf5.html>. Both snapshot formats contain exactly the same information. I.e. the conversion is lossless.\n",
    "\n",
    "Plonk aims to be neutral to the SPH code that generated the data. It is possible to add readers for other SPH codes. However, if the file format is non-standard, i.e. not HDF5 or similar, it is non-trivial to add support."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting started\n",
    "\n",
    "I assume you have Plonk installed and in a working state. If not see the documentation. If you are already using conda-forge (`conda config --add channels conda-forge`), it should be as easy as \n",
    "\n",
    "```\n",
    "conda install plonk\n",
    "```\n",
    "\n",
    "To get started we will download a data set from a Phantom simulation of a dust and gas protoplanetary disc with an embedded protoplanet, and perform some basic tasks.\n",
    "\n",
    "### Download sample data\n",
    "\n",
    "Download the sample data set, `plonk_example_data.tar`, from [Anaconda Cloud]( https://anaconda.org/dmentipl/plonk_example_data/files). (Note: the data set is approximately 850 MB.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget \"https://anaconda.org/dmentipl/plonk_example_data/1/download/plonk_example_data.tar\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then extract the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tar xvf plonk_example_data.tar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change into the `plonk_example_data` directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd plonk_example_data/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with SPH data\n",
    "\n",
    "Here we demonstrate some basic usage of Plonk to read a Phantom SPH data set, get some particle arrays as NumPy arrays, and perform some basic analysis and visualization.\n",
    "\n",
    "First import Plonk, matplotlib, and NumPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plonk\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can check the version of Plonk to make sure you're up to date. (It should be `0.3.0`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plonk.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take advantage of the help function. All functions, classes, modules, and sub-packages in Plonk have built-in documentation, including, at least, instructions on how to call the function or instantiate the object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(plonk.Snap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Snapshots\n",
    "\n",
    "SPH snapshot files are represented by the `Snap` class. This object contains a properties dictionary, particle arrays, which are lazily loaded from file. Here we demonstrate instantiating a `Snap` object, and accessing some properties and particle arrays.\n",
    "\n",
    "First, we load the snapshot with the `load_snap()` function. You can pass a string or `pathlib.Path` object to point to the location of the snapshot in the file system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'disc_00030.h5'\n",
    "snap = plonk.load_snap(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can access arrays by their name passed in as a string."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snap['position']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There may be a small delay as the data is read from file. After the array is read from file it is cached in memory, so that subsequent calls are faster.\n",
    "\n",
    "To see what arrays are loaded into memory you can use the `loaded_arrays()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snap.loaded_arrays()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `available_arrays()` to see what arrays are available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snap.available_arrays()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also define your own alias to access arrays. For example, if you prefer to use the name ‘coordinate’ rather than ‘position’, use the `add_alias()` method to add an alias."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snap.add_alias(name='position', alias='coordinate')\n",
    "snap['coordinate']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Snap` object has a `properties` attribute which is a dictionary of metadata, i.e. non-array data, on the snapshot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snap.properties['time']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(snap.properties)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sink particles are handled separately from the fluid, e.g. gas or dust, particles. They are available as an attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snap.sinks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snap.sinks['spin']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation\n",
    "\n",
    "SPH simulation data is usually spread over multiple files of, possibly, different types, even though, logically, a simulation is a singular “object”. Plonk has the `Simulation` class to represent the complete data set. `Simulation` is an aggregation of the `Snap` and `Evolution` (see below) objects, plus metadata, such as the directory on the file system.\n",
    "\n",
    "Use the `load_sim()` function to instantiate a `Simulation` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prefix = 'disc'\n",
    "sim = plonk.load_sim(prefix=prefix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of the snapshots are available via `snaps` as a list. We can get the first five snapshots with the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.snaps[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Simulation` class has attributes `global_quantities` and `sink_quantities` which are instances of the `Evolution` class, discussed in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evolution\n",
    "\n",
    "SPH simulation data also include auxiliary files containing globally-averaged quantities output more frequently than snapshot files. For example, Phantom writes text files with the suffix `.ev`. These files are output every time step rather than at the frequency of the snapshot files.\n",
    "\n",
    "The Plonk `Evolution` class encapsulates this data. Use `load_ev()` to instantiate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev = plonk.load_ev('disc01.ev')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data may be split over several files, for example, if the simulation was run with multiple jobs on a computation cluster. In that case, pass in a tuple or list of files in chronological order to `load_ev()`, and Plonk will concatenate the data removing any duplicated time steps.\n",
    "\n",
    "The underlying data is stored as a pandas DataFrame. This allows for the use of typical pandas operations with which users in the scientific Python community may be familiar with. Access the DataFrame with the `data` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can plot columns with the pandas plotting interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev.plot('time', ['xcom', 'ycom', 'zcom'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization of SPH data\n",
    "\n",
    "SPH particle data is not gridded like the data produced by, for example, finite difference or finite volume hydrodynamical codes. One visualization method is to plot the particles as a scatter plot, and possibly color the particles with the magnitude of a quantity of interest. An alternative is to interpolate any quantity on the particles to a pixel grid with weighted kernel density estimation. For the technical details, see [Price (2007), PASA, 24, 3, 159](https://ui.adsabs.harvard.edu/abs/2007PASA...24..159P).\n",
    "\n",
    "You can use the `visualize.render()` function to interpolate a quantity to a pixel grid to show as an image. For example, in the following we produce a plot of column density, i.e. a projection plot.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viz = plonk.visualize.render(snap=snap, quantity='density', extent=[-200, 200, -200, 200])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This produces an image via Matplotlib. The `Visualization` object has attributes to access the underlying Matplotlib objects, and their methods. For example, we can use the Matplotlib `image.AxesImage` object to set the limits of the colorbar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viz.image.set_clim(vmin=0.5e-8, vmax=1.5e-8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NB: the previous line won't do anything if the image in displayed inline in the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extra options can be passed in via dictionaries. For example, we can produce a polar plot with the scalar_options keyword argument as follows. We also change the colormap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viz = plonk.visualize.render(\n",
    "    sim.snaps[-5],\n",
    "    quantity='rho',\n",
    "    scalar_options={'polar_coordinates': True, 'cmap': 'viridis'},\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More fine-grained control can be achieved by using the function `visualize.plot()`. See the [API documentation](https://plonk.readthedocs.io/en/latest/api.html) for more details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of SPH data\n",
    "\n",
    "### Subsnaps\n",
    "\n",
    "When analyzing SPH data it can be useful to look at a subset of particles. For example, the simulation we have been working with has dust and gas. So far we have been plotting the total density. We may want to visualize the dust and gas separately.\n",
    "\n",
    "To do this we take a `SubSnap`. The ‘dust_id’ array distinguishes between dust and gas particles. Gas particles have a ‘dust_id’ of 0. Dust particles have a ‘dust_id’ of 1 (or greater for multiple species). In this simulation there is only one dust species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gas = snap[snap['dust_id'] == 0]\n",
    "dust = snap[snap['dust_id'] == 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can access arrays on the `SubSnap` objects as for any `Snap` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gas['mass']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dust['mass']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s plot the gas and dust side-by-side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subsnaps = [gas, dust]\n",
    "extent = (-200, 200, -200, 200)\n",
    "\n",
    "fig, axes = plt.subplots(ncols=2, figsize=(12, 5))\n",
    "\n",
    "for subsnap, axis in zip(subsnaps, axes):\n",
    "    plonk.visualize.render(\n",
    "        subsnap, quantity='rho', extent=extent, axis=axis\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derived arrays\n",
    "\n",
    "Sometimes you need new arrays on the particles that are not available in the snapshot files. You can create a new, derived array on the particles as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snap['r'] = np.sqrt(snap['x'] ** 2 + snap['y'] ** 2)\n",
    "snap['r']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where, here, we have used the fact that Plonk knows that ‘x’ and ‘y’ refer to the x- and y-components of the position array.\n",
    "\n",
    "Alternatively, you can define a function for a derived array. This makes use of the decorator `add_array()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@plonk.Snap.add_array\n",
    "def radius(snap):\n",
    "    radius = np.hypot(snap['x'], snap['y'])\n",
    "    return radius"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snap['radius']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Profiles\n",
    "\n",
    "Generating a radial profile is a convenient method to reduce the dimensionality of the full data set. For example, we may want to see how the surface density and aspect ratio of the disc vary with radius.\n",
    "\n",
    "To do this we use the `Profile` class in the `analysis` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prof = plonk.analysis.Profile(snap, radius_min=10, radius_max=200)\n",
    "prof"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see what profiles are loaded and what are available use the `loaded_keys()` and `available_keys()` methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prof.loaded_keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prof.available_keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To load a profile, simply call it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prof['scale_height']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can convert the data in the `Profile` object to a pandas DataFrame with the `to_dataframe()` method. This takes all loaded profiles and puts them into the DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "profiles = (\n",
    "   'angmom_mag',\n",
    "   'angmom_phi',\n",
    "   'angmom_theta',\n",
    "   'density',\n",
    "   'scale_height',\n",
    ")\n",
    "for p in profiles:\n",
    "    prof[p]\n",
    "df = prof.to_dataframe()\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can use pandas plotting methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('seaborn'):\n",
    "    fig, axes = plt.subplots(ncols=2, figsize=(12, 5))\n",
    "    df.plot('radius', 'density', ax=axes[0])\n",
    "    df.plot('radius', 'scale_height', ax=axes[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrapping up\n",
    "\n",
    "- To see more examples, go to <https://plonk.readthedocs.io/en/latest/examples.html>.\n",
    "- If you find a bug, please raise an issue on GitHub at <https://github.com/dmentipl/plonk/issues>.\n",
    "- If Plonk is missing a feature you want, you can raise a \"feature request\" issue on GitHub.\n",
    "- Alternatively, you can contribute your own code to Plonk, see [CONTRIBUTING.md](https://github.com/dmentipl/plonk/blob/master/CONTRIBUTING.md) for details. Contributions are welcome 😃"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}