04 - Running_an_Nbody_code.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-ph4 amuse-huayno\"\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": {
    "tags": []
   },
   "source": [
    "How to run a simple N-body code\n",
    "====================\n",
    "\n",
    "Here we will generate initial conditions for an N-body code, run a small simulation and analyse the results. This analysis is performed on a 100 star cluster in a 1 pc virial-radius King model. Stellar masses are taken randomly from a Salpeter distribution.\n",
    "Stellar evolution is ignored in this simulation.\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 initial conditions (mass function, King-model)\n",
    "- to initialize a gravitational N-body code\n",
    "- channels and intra-code data transfer\n",
    "- detecting binaries\n",
    "- plotting results\n",
    "- making cumulative distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy\n",
    "from matplotlib import pyplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.units import units, constants\n",
    "from amuse.lab import new_powerlaw_mass_distribution\n",
    "\n",
    "from amuse.units import nbody_system\n",
    "from amuse.ic.kingmodel import new_king_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by setting-up a simulation by specifying a stellar mass distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_stars = 100\n",
    "alpha_IMF = -2.35\n",
    "\n",
    "m_stars = new_powerlaw_mass_distribution(\n",
    "    n_stars, 0.1 | units.MSun, \n",
    "    100 | units.MSun, alpha_IMF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now check to see if this indeed gives one the expected mass function by plotting the cumulative distribution against a theoretical distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_cdf(\n",
    "    m, alpha_IMF=-2.35):\n",
    "    m=m.sorted()\n",
    "    fm=m.cumsum()\n",
    "    logm_max=numpy.log10(m[-1].value_in(units.MSun))\n",
    "    x = 10**numpy.arange(-1.0, logm_max, 1/100)\n",
    "    y = x**(alpha_IMF+2)\n",
    "    fy = numpy.cumsum(y)\n",
    "    pyplot.semilogx(x, fy/fy[-1], lw=2, ls=\":\")\n",
    "    pyplot.semilogx(m.value_in(units.MSun), fm/fm[-1], lw=4, ls=\"-\")\n",
    "\n",
    "plot_cdf(m_stars, alpha_IMF=alpha_IMF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question:\n",
    "Can you explain why the dotted curve does not overplot with the solid curve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can initialize the cluster size and declare the converter.\n",
    "\n",
    "### Converter\n",
    "N-body codes often operate using special dimensionless N-body units (where the gravitational constant is set to unity). This is convenient because it means that the results can be scaled to different systems. AMUSE can use the dimensionless unit system of the codes, using the module nbody_system (which defines the necessary unit system)\n",
    "\n",
    "However, practically speaking as soon as we have a physical scale in the system, such as here, the minimal and maximum mass of the IMF we want our simulation to scale to a definite unit system, and it is better to simply use SI units. In order to make the code understand SI units, we must pass it a unit converter. This contains two quantities in the relevant unit system, from which AMUSE derives the proper unit conversion. Best practice is to keep those quantities on the order of the system's scale (e.g. the total mass and the radius of the outermost planet in a planetary system) to ensure numerical stability.\n",
    "\n",
    "Also see the notebook tutorial on generic units on the importance of converters in AMUSE. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_cluster = 1.0 | units.parsec\n",
    "converter = nbody_system.nbody_to_si(m_stars.sum(), r_cluster)\n",
    "\n",
    "W0 = 3.0\n",
    "bodies = new_king_model(n_stars, W0, convert_nbody=converter)\n",
    "bodies.move_to_center()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned above, the nbody system is defined by having the gravitational constant G have the numerical value 1. This leaves the actual unit system for purposes of dynamical system calculations underdetermined (because the unit space is spanned by three dimensions, e.g. time, mass and length). A *converter* provides a way to scale or convert units. A converter is instantiated with two independent physical quantities (which must also be independent of G) that are assumed to scale such that the numerical value is 1 in the corresponding nbody_system dimensions. Together with G=1 this defines a mapping between nbody_system and SI. It is also often times passed to some initial condition generator, in which case it then implicitly scales the system!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"for this particular converter:\")\n",
    "print(\n",
    "    f\"  The gravitational constant G in nbody units is: {converter.to_nbody(constants.G)}. \"\n",
    "    \"Note that this is not unitless!\")\n",
    "print(f\"  1 Myr is {converter.to_nbody(1. | units.Myr)} in nbody units\")\n",
    "print(f\"  1 nbody time corresponds to {converter.to_si(1 | nbody_system.time)} in SI units\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the result by plotting the X-Y positions of the bodies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.plot import scatter\n",
    "\n",
    "v = (bodies.vx**2 + bodies.vy**2 + bodies.vz**2).sqrt()\n",
    "scatter(bodies.x, bodies.y, c=v.value_in(units.kms), alpha=0.5)\n",
    "pyplot.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code evolves the cluster in time by using one of AMUSEs nbody integrators. We will explain a bit more below, but try to read see if you can understand the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.community.ph4.interface import Ph4\n",
    "from amuse.ext.LagrangianRadii import LagrangianRadii\n",
    "\n",
    "gravity = Ph4(converter)\n",
    "gravity.particles.add_particles(bodies)\n",
    "channel = gravity.particles.new_channel_to(bodies)\n",
    "\n",
    "\n",
    "times = numpy.arange(0, 100, 0.1) | units.Myr\n",
    "RL25 = [] | units.parsec\n",
    "Rvir = [] | units.parsec\n",
    "for time in times:\n",
    "    gravity.evolve_model(time)\n",
    "    channel.copy()  # Copy from gravity.particles to bodies\n",
    "    Rvir.append(bodies.virial_radius())\n",
    "    L = LagrangianRadii(bodies)\n",
    "    RL25.append(LagrangianRadii(bodies)[5])\n",
    "\n",
    "    if not time.value_in(units.Myr) % 10.0:\n",
    "        print(\"cluster at Time=\", time.in_(units.Myr), \n",
    "              \"Mass=\", bodies.mass.sum().in_(units.MSun),\n",
    "              \"Rvir=\", Rvir[-1].in_(units.parsec))\n",
    "    b = bodies.get_binaries()\n",
    "    if(len(b) > 0):\n",
    "        print(\"Number of binaries found:\", len(b))\n",
    "        print(\"first binary:\", b[0])\n",
    "\n",
    "gravity.stop()\n",
    "        \n",
    "pyplot.plot(times.value_in(units.Myr), RL25.value_in(units.parsec))\n",
    "pyplot.plot(times.value_in(units.Myr), Rvir.value_in(units.parsec))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have used a channel. A channel is a 'permanent' connection between two particle sets. In this case, between the particle set in the code and the original\n",
    "one, generated by the initial condition procedure. Updates are effected by the copy() method - this copies over all attributes of the set (overwriting the existing data).\n",
    "So each copy updates the particle positions and velocities in bodies with evolution as calculated by the code.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have now calculated the dynamical evolution of a small cluster of stars including a realistic mass function.\n",
    "\n",
    "Assignments and questions:\n",
    "---------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 1\n",
    "\n",
    "The AMUSE interface for N-body code is defined [here](https://github.com/amusecode/amuse/blob/main/src/amuse/community/interface/gd.py]). While it is not within the scope of this tutorial to fully learn this. This would also be a very unnecessary because most of the interface functions are hidden! A few elements are important: try out and examine the following properties and methods (what is the difference between methods and properties?) :\n",
    "`kinetic_energy`, `potential_energy`, `total_energy`, `get_potential_at_point`, `get_gravity_at_point`, `model_time`\n",
    "and examine the `parameters` and `particles` attributes. Try to see how they behave if you add particles and then remove again etc.\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.ic.plummer import new_plummer_model\n",
    "from amuse.community.huayno.interface import Huayno\n",
    "\n",
    "p = new_plummer_model(100)\n",
    "h = Huayno()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 1:\n",
    "Run the notebook and describe the evolution of the virial radius."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 2:\n",
    "There is a check for binary detection in the script.\n",
    "If no binary formed in your run, that is okay.\n",
    "\n",
    "Run the script again until a binary forms. Note that bodies is updated during the run! If you want to speed up binary formation, you can remake the initial conditions and\n",
    "slow down all the stars by multiplying the velocity by something <1, say 0.9. Observe the cange in the evolution of the virial radius in this case!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 2:\n",
    "One of the parameters in the binary is the *hardness*.\n",
    "What is the hardness of the binary from your last run, and explain what is meant with this parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 3:\n",
    "Rewrite the script in such a way that the event loop over time continues indefinitely, but that the code stops as soon as the first binary is detected.\n",
    "\n",
    "Run the script several times until you have enough first binaries formed to make a histogram of the first-binary formation time-scale.\n",
    "\n",
    "Now, make a small change to your initial conditions.\n",
    "A variation could include the dimensionless depth of the potential of the initial King model, the slope of the initial-mass function, its lower limit, or the number of stars in the simulation).\n",
    "\n",
    "Now redo the calculation in which you generate a histogram of moments of first-binary formation.\n",
    "\n",
    "Make one figure with both distributions: the original distribution for first-binary formation time and the second series of runs (with the altered initial conditions.\n",
    "Present both distributions as histograms and as cumulative distributions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 3:\n",
    "Can you understand the difference in first-binary formation-time based on the changes you introduced in the initial conditions?\n",
    "\n",
    "Did you perform enough runs to make a statistically significant statement about the time of first-binary formation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 4:\n",
    "Make a cumulative distribution of the masses of the two stars for each first binary and compare it to the initial mass-function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 4:\n",
    "Can you understand the difference between the typical masses of the binaries that formed first and the initial mass-function?"
   ]
  },
  {
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}