07 - Running_Nbody_with_collisions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "# preamble script to check and install AMUSE components if necessary\n",
    "\n",
    "# required packages for this tutorial:\n",
    "PACKAGES=\"mpi4py amuse-framework amuse-seba amuse-sse amuse-ph4\"\n",
    "# skip in case a full development install is present\n",
    "pip show amuse-devel && exit 0\n",
    "for package in ${PACKAGES} \n",
    "do\n",
    "  pip show ${package} || pip install ${package}\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# the following fixes are highly recommended\n",
    "\n",
    "#allow oversubscription for openMPI\n",
    "import os\n",
    "os.environ[\"OMPI_MCA_rmaps_base_oversubscribe\"]=\"true\"\n",
    "\n",
    "# use lower cpu resources for idle codes\n",
    "from amuse.support import options\n",
    "options.GlobalOptions.instance().override_value_for_option(\"polling_interval_in_milliseconds\", 10)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How to run a code with stopping conditions\n",
    "====================\n",
    "\n",
    "Stopping conditions should be considered interrupts that are initiated through special conditions in the running code. Stopping conditions leave the running code in a restartable state while exiting the code to hand the control back to the python (i.e. the AMUSE) script. There are stopping conditions in most codes, but here we only discuss one of them in the stellar dynamics codes, which is called: *collision_detection*\n",
    "\n",
    "For reference you can read chapter 2 of Portegies Zwart & McMillan 2018 (2018araa.book.....P). \n",
    "\n",
    "With this tutorial you will learn\n",
    "- how to generate inital conditions (mass-function, King-model)\n",
    "- to initialize an gravitational N-body code\n",
    "- channels and intra-code data transfer\n",
    "- detecting collisions\n",
    "- plotting results\n",
    "- making Hertzsprung-Russell diagrams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy\n",
    "numpy.random.seed(98429435)\n",
    "from matplotlib import pyplot\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by generating the initial conditions and initializing the N-body code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from amuse.units import units\n",
    "from amuse.units import nbody_system\n",
    "from amuse.ic.plummer import new_plummer_model\n",
    "from amuse.lab import new_powerlaw_mass_distribution\n",
    "from amuse.lab import Particles\n",
    "\n",
    "from amuse.community.ph4 import Ph4\n",
    "from amuse.community.seba import Seba\n",
    "from amuse.ext.LagrangianRadii import LagrangianRadii"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n_stars = 100\n",
    "alpha_IMF = -2.35\n",
    "\n",
    "m_stars = new_powerlaw_mass_distribution(\n",
    "    n_stars,\n",
    "    10.0 | units.MSun,\n",
    "    100.0 | units.MSun,\n",
    "    alpha_IMF)\n",
    "\n",
    "r_cluster = 1.0 | units.parsec\n",
    "converter = nbody_system.nbody_to_si(m_stars.sum(), r_cluster)\n",
    "\n",
    "stars = new_plummer_model(n_stars, convert_nbody=converter)\n",
    "stars.mass = m_stars\n",
    "stars.collision_radius= 0. | units.RSun \n",
    "stars.scale_to_standard(converter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "stellar = Seba()\n",
    "stellar.particles.add_particles(stars)\n",
    "\n",
    "gravity = Ph4(converter, number_of_workers=6)\n",
    "gravity.particles.add_particles(stars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's declare the channels.\n",
    "This is somewhat more elaborate here, because there is a new parameter to declare, which is *collision_radius*.\n",
    "The parameter is needed in order to ensure that our stellar collisions behave as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "stellar_attributes = [\"mass\", \"radius\", \"age\", \"temperature\", \"luminosity\"]\n",
    "channel = {\"from_stellar\":\n",
    "              stellar.particles.new_channel_to(stars,\n",
    "              attributes=stellar_attributes,\n",
    "              target_names=stellar_attributes),\n",
    "            \"from_gravity\": \n",
    "              gravity.particles.new_channel_to(stars,\n",
    "              attributes=[\"x\", \"y\", \"z\", \"vx\", \"vy\", \"vz\", \"mass\"],\n",
    "              target_names=[\"x\", \"y\", \"z\", \"vx\", \"vy\", \"vz\", \"mass\"]),\n",
    "            \"to_gravity\": \n",
    "              stars.new_channel_to(gravity.particles,\n",
    "              attributes=[\"mass\", \"collision_radius\"],\n",
    "              target_names=[\"mass\", \"radius\"])}\n",
    "    \n",
    "channel[\"from_stellar\"].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We here explicitely copied particle attributes which we think are of interesting for this experiment, including stellar mass, radius, age, temperature and luminosity.\n",
    "\n",
    "##### *collision_radius*\n",
    "An additional parameter was introduced to keep track of the collisional radius for the stars. Stars tend to be small compared to the inter-stellar distance. In order to warrant some nice collisions in our simulation we may want to expand the stellar size in comparison with the equilibrium photospheric radius of the stars. We can do this by adding an additional parameter, here called *collision_radius* to the particle set. Notice that the *collision_radius* is copied into the attribute *radius* of the stars in the N-body code, because there is no separate collision-radius in the N-body code, but the attribute *radius* is used to check for two objects to come close before an exception is thrown."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def plot_snapshot(bodies):\n",
    "    v = (bodies.vx**2 + bodies.vy**2 + bodies.vz**2).sqrt()\n",
    "    s = bodies.mass.value_in(units.MSun)\n",
    "    fig, (ax1, ax2) = pyplot.subplots(1, 2, figsize=(11, 3))\n",
    "    ax1.scatter(stars.temperature.value_in(units.K),\n",
    "                stars.luminosity.value_in(units.LSun), \n",
    "                c=v.value_in(units.kms),\n",
    "                s=s)\n",
    "    ax1.set_xlim(6.e+4, 20000)\n",
    "    ax1.set_ylim(1.e+3, 1.e+7)\n",
    "    ax1.loglog()\n",
    "    ax1.set_xlabel(\"T [K]\")\n",
    "    ax1.set_ylabel(\"L [$L_\\odot$]\")\n",
    "    ax2.scatter(bodies.x.value_in(units.pc), \n",
    "                bodies.y.value_in(units.pc), \n",
    "                c=v.value_in(units.kms), \n",
    "                s=s)\n",
    "    pyplot.gca().set_aspect('equal', adjustable='box')\n",
    "    ax2.set_xlabel(\"x [pc]\")\n",
    "    ax2.set_ylabel(\"y [pc]\")\n",
    "    ax2.set_xlim(-5, 5)\n",
    "    ax2.set_ylim(-5, 5)\n",
    "    pyplot.show()\n",
    "plot_snapshot(stars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now initiate the stopping conditions in order to detect a merger between two stars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "stopping_condition = gravity.stopping_conditions.collision_detection\n",
    "stopping_condition.enable()\n",
    "collision_radius_multiplication_factor = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We introduced a multiplication factor in order to acquire many nice collisions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def merge_two_stars(bodies, particles_in_encounter):\n",
    "    com_pos = particles_in_encounter.center_of_mass()\n",
    "    com_vel = particles_in_encounter.center_of_mass_velocity()\n",
    "    d = (particles_in_encounter[0].position - particles_in_encounter[1].position)\n",
    "    v = (particles_in_encounter[0].velocity - particles_in_encounter[1].velocity)\n",
    "    print(\"Actually merger occurred:\")\n",
    "    print(\"Two stars (M=\",particles_in_encounter.mass.in_(units.MSun),\n",
    "          \") collided with d=\", d.length().in_(units.au))\n",
    "    new_particle=Particles(1)\n",
    "    new_particle.mass = particles_in_encounter.total_mass()\n",
    "    new_particle.age = min(particles_in_encounter.age) \\\n",
    "                     * max(particles_in_encounter.mass)/new_particle.mass\n",
    "    new_particle.position = com_pos\n",
    "    new_particle.velocity = com_vel\n",
    "    new_particle.radius = particles_in_encounter.radius.sum()\n",
    "    bodies.add_particles(new_particle)\n",
    "    bodies.remove_particles(particles_in_encounter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When two stars collide, we will merge with the function above. We create a new particle with the summed mass of the two merged stars, and remove the old stars. What quantities are conserved during this merger process?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def resolve_collision(collision_detection, gravity, stellar, bodies):\n",
    "    if collision_detection.is_set():\n",
    "        print(\"Well, we have an actual collision between two or more stars.\")\n",
    "        print(\"This happened at time=\", stellar.model_time.in_(units.Myr))\n",
    "        for ci in range(len(collision_detection.particles(0))): \n",
    "            encountering_particles = Particles(particles=[collision_detection.particles(0)[ci],\n",
    "                                                          collision_detection.particles(1)[ci]])\n",
    "            colliding_stars = encountering_particles.get_intersecting_subset_in(bodies)\n",
    "            merge_two_stars(bodies, colliding_stars)\n",
    "            bodies.synchronize_to(gravity.particles)\n",
    "            bodies.synchronize_to(stellar.particles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The collision detector is a separate object from the code. If this is triggered, we process all pairs of stars that have collided."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "end_time = 10.0 | units.Myr\n",
    "model_time = 0 | units.Myr\n",
    "rvir = [] | units.pc\n",
    "time = [] | units.Myr\n",
    "t_diag = 1| units.Myr\n",
    "while(model_time<end_time):\n",
    "    dt = stellar.particles.time_step.min()\n",
    "    model_time += dt\n",
    "    stellar.evolve_model(model_time)\n",
    "    channel[\"from_stellar\"].copy()\n",
    "    stars.collision_radius = stars.radius * collision_radius_multiplication_factor\n",
    "    channel[\"to_gravity\"].copy()\n",
    "    gravity.evolve_model(model_time)\n",
    "    resolve_collision(stopping_condition, gravity, stellar, stars)\n",
    "    channel[\"from_gravity\"].copy()\n",
    "    \n",
    "    if model_time >= t_diag:\n",
    "        t_diag += 1 | units.Myr\n",
    "\n",
    "        print(\"Evolved to t=\", stellar.model_time.in_(units.Myr), \n",
    "              gravity.model_time.in_(units.Myr),\n",
    "              \"N=\", len(stars),\n",
    "              \"mass=\", stars.mass.sum().in_(units.MSun), \n",
    "              \"rvir=\", stars.virial_radius().in_(units.pc))\n",
    "    rvir.append(stars.virial_radius())\n",
    "    time.append(model_time)\n",
    "    \n",
    "print(\"done\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "stellar.stop()\n",
    "gravity.stop()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_snapshot(stars)\n",
    "pyplot.plot(time.value_in(units.Myr), rvir.value_in(units.parsec))\n",
    "pyplot.xlabel(\"t [Myr]\")\n",
    "pyplot.ylabel(\"$R_{vir}$ [pc]\")\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have calculated the evolution of a small cluster of 100 stars. Your calculation included Newtonian dynamics, stellar evolution and collisions.\n",
    "The collisions were rather simple, sticky spheres, and we have chosen the stellar radii to be a 1000 times larger than normal.\n",
    "\n",
    "Assignments and questions:\n",
    "---------------\n",
    "\n",
    "### Question 1:\n",
    "\n",
    " * Why did we run with a random seed?\n",
    " * How many collisions did occur?\n",
    " * What are these curious stars that are not on the main sequence in the Hertzsprung-Russel diagram?\n",
    " * Where in the cluster (with respect to the core-radius and virial-radius) are the collision products? Note that the core radius can be acquired at runtime from the particle set via the member function: *stars.densitycentre_coreradius_coredens* (which returns the position vector of the density centre, the core radius, and the core density).\n",
    " \n",
    "### Assignment 1:\n",
    "Include the appropriate diagnostics in the script in order to address the following questions:\n",
    " \n",
    " * What was the distance from the cluster center where the collissions took place?\n",
    " * What is the amount of energy that has been dissipated in the collisions?\n",
    " * What is the actual numerical energy error in the calculation?\n",
    " \n",
    "### Assignment 2:\n",
    "We assumed that mass and momentum are conserved in the collision process. Change the script in such a way that the mass of the collision product is the sum of the primary mass and half that of the secondary.\n",
    "Rerun the script and explain the differences.\n",
    "\n",
    "### Questions 2:\n",
    "\n",
    " * Has the change affected the number of collisions?\n",
    " * What change do you see in the Hertzsprung-Russel diagram?\n",
    " * What happened to the cluster virial radius?\n",
    " \n",
    "### Questions 3:\n",
    "Considerable magic seems to be happening in the routine *resolve_collision*.\n",
    "\n",
    " * explain what happens in the following lines:\n",
    "```python\n",
    "   encountering_particles = Particles(\n",
    "       particles=[collision_detection.particles(0)[ci],\n",
    "       collision_detection.particles(1)[ci]])\n",
    "```\n",
    " \n",
    " * and what happens in:\n",
    "```python\n",
    "colliding_stars = encountering_particles.get_intersecting_subset_in(bodies)\n",
    "```\n",
    "\n",
    " * and in:\n",
    "```python\n",
    "bodies.synchronize_to(gravity.particles)\n",
    "bodies.synchronize_to(stellar.particles)\n",
    "```\n",
    "\n",
    "### Assignment 4:\n",
    " * Write a one-line python script to find the brighest star."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}