running_hydrodynamics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How to run a smoothed particles Hydrodynamics simulation\n",
    "====================\n",
    "\n",
    "We will run a hydrodynamical simulation using a smoothed-particles hydrodynamics code.\n",
    "\n",
    "For reference you can read chapter 5 of Portegies Zwart & McMillan 2018 (2018araa.book.....P). \n",
    "\n",
    "For this experiment, we will evolve a high-mass star to the end of its lifetime using a Henyey stellar evolution code (this may take a while).\n",
    "We subsequently convert the star to a hydrodynamical realization in particles and inject the supernova energy to the inner few SPH-particles. We further simulate the evolution of the gas-blob in time using a hydrodynamics code.\n",
    "\n",
    "To reduce the runtime of the calculation, we do not evolve the star to the end of its evolution, but just for a short while. Of course, the exploding star, in that case, will still be that of a main-sequence star.\n",
    "\n",
    "The script in this notebook is based on the work by [Rimoldi et al.](https://ui.adsabs.harvard.edu/abs/2016ComAC...3....2R/abstract).\n",
    "\n",
    "With this tutorial you will learn\n",
    " * Run another AMUSE module to generate initial conditions for yet another code.\n",
    " * How to recover the crash of a code and pick up the result.\n",
    " * Store simulation data in the form of python pickle files, and recover from those.\n",
    " * Plot the result of a hydrodynamical simulation.\n",
    " * Make an animation of simulation results.\n",
    " * Run an AMUSE module as a parallel job.\n",
    " * Evolve a star to a certain evolutionary stage.\n",
    " \n",
    "This tutorial will require the stellar evolution code MESA to evolve the star. A separate file with a pre-calculated stellar sturcture is available (super_giant_stellar_structure.pkl), but reading in this file still requires MESA to be installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import numpy\n",
    "from matplotlib import pyplot\n",
    "from amuse.units import units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After loading in the most essential modules, we will set up the stellar evolution. In case MESA is unavailable, the repository contains a copy of the resulting stellar structure file. Simply switch which line defining *pickle_file* is commented out. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.ext.star_to_sph import (pickle_stellar_model, convert_stellar_model_to_SPH,)\n",
    "from amuse.test.amusetest import get_path_to_results\n",
    "from amuse.community.mesa.interface import MESA\n",
    "from amuse.datamodel import Particles\n",
    "\n",
    "def setup_stellar_evolution_model():\n",
    "    out_pickle_file = os.path.join(get_path_to_results(), \n",
    "                                   \"super_giant_stellar_structure.pkl\")\n",
    "\n",
    "    stellar_evolution = MESA(redirection=\"none\")\n",
    "    stars = Particles(1)\n",
    "    stars.mass = 15.0 | units.MSun\n",
    "    stellar_evolution.particles.add_particles(stars)\n",
    "    stellar_evolution.commit_particles()\n",
    "\n",
    "    print(\n",
    "        \"Evolving a MESA star with mass:\",\n",
    "        stellar_evolution.particles[0].mass\n",
    "    )\n",
    "    try:\n",
    "        while stellar_evolution.model_time<0.12|units.Myr:\n",
    "            stellar_evolution.evolve_model()\n",
    "            print(\"star:\", stellar_evolution.particles[0].stellar_type, stellar_evolution.model_time.in_(units.Myr))\n",
    "    except AmuseException as ex:\n",
    "        print(\"Evolved star to\", stellar_evolution.particles[0].age)\n",
    "        print(\"Radius:\", stellar_evolution.particles[0].radius)\n",
    "    \n",
    "    pickle_stellar_model(stellar_evolution.particles[0], out_pickle_file)\n",
    "    stellar_evolution.stop()\n",
    "    return out_pickle_file\n",
    "\n",
    "pickle_file = setup_stellar_evolution_model()\n",
    "#pickle_file = './super_giant_stellar_structure.pkl'\n",
    "print(\"Star generated.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now set up the stellar model in a Henyey-code compatible form. Now we have to convert that into a distribution of SPH particles. There is a function in AMUSE to do this. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_sph_particles = 100\n",
    "print(pickle_file)\n",
    "print(\"Creating initial conditions from a MESA stellar evolution model...\")\n",
    "model = convert_stellar_model_to_SPH(\n",
    "        None,\n",
    "        number_of_sph_particles,\n",
    "        seed=12345,\n",
    "        pickle_file=pickle_file,\n",
    "        #        base_grid_options = dict(type = \"glass\", target_rms = 0.01),\n",
    "        with_core_particle=True,\n",
    "        target_core_mass = 1.4|units.MSun\n",
    "    )\n",
    "print(\"model=\", model)\n",
    "core, gas_without_core, core_radius = \\\n",
    "        model.core_particle, model.gas_particles, model.core_radius\n",
    "print(\"Created\", len(gas_without_core),\n",
    "       \"SPH particles and one 'core-particle':\\n\", core)\n",
    "print(\"Setting gravitational smoothing to:\", core_radius.in_(units.km))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_star(core, gas):\n",
    "    pyplot.scatter(gas.x.value_in(units.au), gas.y.value_in(units.au), c='b')\n",
    "    pyplot.scatter(core.x.value_in(units.au), core.y.value_in(units.au), c='r')\n",
    "    pyplot.axis(\"equal\")\n",
    "    pyplot.show()\n",
    "\n",
    "plot_star(core, gas_without_core)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the mass of the core is not precisely what you put in. \n",
    "The reason is that mass has to be preserved in the simulation, and we remove a number of SPH particles to accomodate the mass of the core. It happens to be the case that the sum of the masses of these SPH particles slightly exceeds our intended core mass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inject_supernova_energy(gas_particles, \n",
    "                            explosion_energy=1.0e+51|units.erg,\n",
    "                            exploding_region=10|units.RSun):\n",
    "    inner = gas_particles.select(\n",
    "        lambda pos: pos.length_squared() < exploding_region**2,\n",
    "        [\"position\"])\n",
    "    print(len(inner), \"innermost particles selected.\")\n",
    "    print(\"Adding\", explosion_energy / inner.total_mass(), \"of supernova \" \\\n",
    "        \"(specific internal) energy to each of the n=\", len(inner), \"SPH particles.\")\n",
    "    inner.u += explosion_energy / inner.total_mass()\n",
    "    \n",
    "inject_supernova_energy(gas_without_core, exploding_region=1|units.RSun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we injected 10^51 erg in the inner-most SPH particles (those within 1 RSun of the stellar center). It turns out that we spread the supernova energy among the inner-most SPH particles. Spreading it over too few particles introduces interesting side effects, such as individual SPH particles being launched, as bullets, from the exploding star. Spreading the supernova energy over too many SPH particles, also has interesting consequences: the explosion would them drive the homologeous expansion of the entire star, rather than an explosion from the inside. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.community.fi.interface import Fi\n",
    "from amuse.lab import nbody_system\n",
    "converter = nbody_system.nbody_to_si(10|units.MSun, core_radius)\n",
    "\n",
    "hydro_code = Fi(converter)\n",
    "hydro_code.parameters.epsilon_squared = core_radius**2\n",
    "hydro_code.parameters.n_smooth_tol = 0.01\n",
    "hydro_code.gas_particles.add_particles(gas_without_core)\n",
    "hydro_code.dm_particles.add_particle(core)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We introduced and initialized the hydrodynamics solver, for which we adopted *Fi*. An alternative would be to use *Gadget2* or a grid-based hydro code.\n",
    "Each hydro solver has its advantages and disadvantages. For SPH, the smoothing length and associated tolerance are important parameters to set. Another important parameter is the artivicial viscosity parameter *alpha*. You may want to experiment with those settings, probably breaking the code on the way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.datamodel import Grid, new_regular_grid\n",
    "\n",
    "def hydro_plot(hydro_code, view_size, npixels):\n",
    "    view = [-1, 1, -1, 1] * view_size\n",
    "    shape = (npixels, npixels, 1)\n",
    "    size = npixels**2\n",
    "    axis_lengths = [0.0, 0.0, 0.0] | units.m\n",
    "    axis_lengths[0] = view[1] - view[0]\n",
    "    axis_lengths[1] = view[3] - view[2]\n",
    "    grid = new_regular_grid(shape, axis_lengths)\n",
    "    grid.x += view[0]\n",
    "    grid.y += view[2]\n",
    "    speed = grid.z.reshape(size) * (0 | 1/units.s)\n",
    "    rho, rhovx, rhovy, rhovz, rhoe = hydro_code.get_hydro_state_at_point(\n",
    "            grid.x.reshape(size),\n",
    "            grid.y.reshape(size),\n",
    "            grid.z.reshape(size), speed, speed, speed)\n",
    "\n",
    "    # we have to make some cuts in the parameter space.\n",
    "    min_v = 800.0 | units.km / units.s\n",
    "    max_v = 3000.0 | units.km / units.s\n",
    "    min_rho = 3.0e-9 | units.g / units.cm**3\n",
    "    max_rho = 1.0e-5 | units.g / units.cm**3\n",
    "    min_E = 1.0e11 | units.J / units.kg\n",
    "    max_E = 1.0e13 | units.J / units.kg\n",
    "\n",
    "    v_sqr = (rhovx**2 + rhovy**2 + rhovz**2) / rho**2\n",
    "    E = rhoe / rho\n",
    "    log_v = numpy.log((v_sqr / min_v**2)) / numpy.log((max_v**2 / min_v**2))\n",
    "    log_rho = numpy.log((rho / min_rho)) / numpy.log((max_rho / min_rho))\n",
    "    log_E = numpy.log((E / min_E)) / numpy.log((max_E / min_E))\n",
    "\n",
    "    red = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum(\n",
    "        numpy.zeros_like(rho.number), log_rho)).reshape(shape)\n",
    "    green = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum(\n",
    "        numpy.zeros_like(rho.number), log_v)).reshape(shape)\n",
    "    blue = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum(\n",
    "        numpy.zeros_like(rho.number), log_E)).reshape(shape)\n",
    "    alpha = numpy.minimum(\n",
    "            numpy.ones_like(log_v),\n",
    "            numpy.maximum(\n",
    "                numpy.zeros_like(log_v),\n",
    "                numpy.log((rho / (10*min_rho)))\n",
    "                )\n",
    "            ).reshape(shape)\n",
    "\n",
    "    rgba = numpy.concatenate((red, green, blue, alpha), axis=2)\n",
    "    \n",
    "    \n",
    "    pyplot.figure(figsize=(npixels/100.0, npixels/100.0), dpi=100)\n",
    "    pyplot.imshow(rgba)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are done setting up the experiment, and can start to run production.\n",
    "In this case, we decided to run in small parts of 10 seconds, 1 minute, 1 hour and 1 day, to see what happens to the exploding star."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydro_code.evolve_model(10.0|units.s)\n",
    "print(\"Done running to time=\", hydro_code.model_time.in_(units.s))\n",
    "view_size = 20 | units.RSun \n",
    "npixels = 200\n",
    "hydro_plot(hydro_code, view_size, npixels)\n",
    "print(\"done plotting\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydro_code.evolve_model(1.0|units.minute)\n",
    "print(\"Done running to time=\", hydro_code.model_time.in_(units.minute))\n",
    "hydro_plot(hydro_code, view_size, npixels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"v-core=\", hydro_code.dm_particles.velocity.length().in_(units.kms))\n",
    "print(\"d-core=\", hydro_code.dm_particles.position.length().in_(units.RSun))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydro_code.evolve_model(1.0 | units.hour)\n",
    "print(\"Done running:\", hydro_code.model_time.in_(units.hour))\n",
    "hydro_plot(hydro_code, view_size, npixels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view_size = 100 | units.RSun\n",
    "while hydro_code.model_time<0.3|units.day:\n",
    "    hydro_code.evolve_model(hydro_code.model_time+(1|units.hour))\n",
    "    print(\"Done running, until:\", hydro_code.model_time.in_(units.hour))\n",
    "    hydro_plot(hydro_code, view_size, npixels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By this time the entire star has probably disappeard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"core=\", hydro_code.dm_particles.velocity.length().in_(units.kms))\n",
    "print(\"core=\", hydro_code.dm_particles.position.length().in_(units.RSun))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hydro_code.stop()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have performed a hydrodynamical simulation of an exploding star. This was realized by first evolving the star to an advanced state using a 1D stellar evolution (a so-called Henyey) code. This one-dimensional structure was subsequently converted to a 3-dimensional hydrodynamical distribution of particles, representing the core and the gaseous envelope.\n",
    "\n",
    "We injected a supernova-amount of energy into the hydrodynamical blob, heating the inner most particles to extreemely high temperatures and let the hydrodynamical code calculate the response.\n",
    "The result was calculated 10 seconds, 1 minute, 1 hour and 1 day after the explosion.\n",
    "\n",
    "\n",
    "Assignments and questions:\n",
    "---------------\n",
    "\n",
    "Note: the script is configured for running with 100 SPH particles, which is too small to study the morphology of the final simulation results, 1 day after the initiation of the supernova. Try running with 1000 or 10000 SPH particles, of better even would be to run with 30000 particles. However, the runs will become rather expensive. Most questions you can answer with ony 1000 SPH particles.\n",
    "\n",
    "### Assignment 1:\n",
    "You can increase the performance of the code by allowing it to run on multiple cores.\n",
    "increase the number of SPH particles from 100 to 1000, and run the code again on as many cores as your computer can muster. Does this indeed increase the code performance? And by how much?\n",
    "\n",
    "Note that Fi has to be set up differently to use multiple cores. Instead of passing *number_of_workers=4* (for 4 workers) when initializing Fi, you have to pass the *mode='openmp'* keyword and set the *OMP_NUM_THREADS* environment variable to the number of workers (using e.g. export or setenv depending on your shell), which has to be done in the same terminal session as your notebook runs in (so you may have to restart). You can check if this worked by also passing Fi the *redirection='none'* keyword. In the terminal running the notebook you should then see *omp threads* and some number.\n",
    "\n",
    "We evolved the star for only 0.12 Myr, which, of course, is too short for a 15 MSun star to even leave the main sequence. A better time scale for such a star is 12Myr.\n",
    "\n",
    "### Assignment 2:\n",
    "Evolve the 15MSun star to 12Myr, and use this stellar structure as input for the hydrodynamics solver.\n",
    "\n",
    "Rather than guessing the main-sequence lifetime of the 15MSun star, one could run the stellar evolution code until it crashes. In that case, the stellar evolution code will automatically terminate while saving the last consistent stellar evolution model. This model can subsequently be used to start the hydrodynamical simulation. This may be preferrable, except that running the stellar evolution code to this point is expensive.\n",
    "\n",
    "In particular the last few steps are time consuming and do not really change that part of the star in which we are most interested (the density gradient between the core-region and the envelope. This part is already converged to its final density profile once the star ascends the second asymptotic giant branch. You can shorten runtime considerably by allowing the star to evolve to that point in its evolution.\n",
    "\n",
    "### Question 1:\n",
    "Where is the stellar model saved, once produced?\n",
    "What is the age and radius of the star at the moment it starts to ascend the first giant branch?\n",
    "\n",
    "### Assignment 3:\n",
    "Rather than evolving the star until the desired moment, we can reuse the same star over and over again without rerunning, at least so long as we do not change the mass of the star. \n",
    "You can use the internal amuse routine *get_path_to_results* to store the picke file in your local amuse repository in $AMUSE_DIR/test_results.\n",
    "\n",
    "Add the necessary lines to the *setup_stellar_evolution_model* script in which the star is evolved to picke the result of the stellar evolution calculation if a new calculation is performed. But first check if the desired pickled star already exist, in which case you should not redo the simulation. \n",
    "\n",
    "In addition, change the script in which the star is evolved to let it run until the star has reached the asymptotic giant branch.\n",
    "\n",
    "### Question 2:\n",
    "In the routine *inject_supernova_energy* we use a *lambda* function. Explain what this function does.\n",
    "\n",
    "### Question 3:\n",
    "In the three images, after running for 10 seconds, 1 minute and 1 hour, the star only seem to change color in the images. Explain why this is.\n",
    "\n",
    "### Assigment 4:\n",
    "In the routine *plot_hydro* we assign various parameters to the colors red, green and blue, and to the transparancy alpha. Explain what these parameters mean, and how this translates to the various colors you observed in the plotted pictures.\n",
    "\n",
    "### Assignment 5:\n",
    "Rerun the calculation but plot the kinetic energy, the potential energy and the total energy of the gas as a function of time from the start of the supernova to 2 days after the supernova. Describe the curves and explain what is happening.\n",
    "\n",
    "### Questions 4:\n",
    "Check the velocity of the neutron star at the end of the simulation.\n",
    "It is not zero, as it was initially. Why is this the case?\n",
    "\n",
    "### Assignment 6:\n",
    "Usually neutron stars are believed to recieve rather high \"natal\" kicks upon birth. We still do not really understand why. In our simulations, we can introduce such a kick velocity by hand. A reasonable kick for a neutron star could be about 250km/s.\n",
    "Impose such a high velocity kick to the newly born neutron star at the moment of the actual explosion and rerung the simulation.\n",
    "\n",
    "### Question 5:\n",
    "What effect that the high natal velocity kick of the neutron star has on the morphology and energy budget of the hydrodynamical outflow of the supernova shell?\n",
    "\n",
    "### Assignment 5:\n",
    "Add the appropriate lines to the plotting routine to write the images to files.\n",
    "Rerun the code, and make an animation of the exploding star."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}