cluster_exercise.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13350a4e-eddd-4f02-bf20-596dc7f5460b",
   "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,
   "id": "c5dcdb48-62b8-4647-b64e-08c15ad9b993",
   "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",
   "id": "37b0881d-c1eb-41fd-b507-f2bd07e96e52",
   "metadata": {
    "tags": []
   },
   "source": [
    "This exercise was originally developed together with Selma de Mink for the UVA computational astrophysics course."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c471a5c",
   "metadata": {},
   "source": [
    "Testing the Evolution of Starclusters with Amuse and making predictions for Gaia \n",
    "================================================================================\n",
    "\n",
    "In this exercise you will be guided to use AMUSE to conduct a numerical n-body experiment and make predictions for Gaia. \n",
    "\n",
    "The bigger question you will address concerns the question of how stars form. We know that stars typically (and possibly always) form in clusters. However, the birth conditions are very hard to measure.  For example: How \"clumpy\" or \"smooth\" are these clusters? How \"virialized\" are they?  \n",
    "\n",
    "You will simulate the evolution of star clusters, with different initial conditions. The aim is to investigate if the initial conditions leave imprints that Gaia can observe.  Specifically you will investigate the population of \"escapers\" and possible \"runaway stars\", i.e. stars that escape their birth cluster.  To do this we will very closely follow the work by Allison (2012, MNRAS 421, 3338). \n",
    "\n",
    "You will report the results of your numerical experiments in the style of a \"letter\" (= short paper normally on a very timely topic.) Detailed instructions below. The questions in this tutorial are to guide you through the experiment. You can use them for the narrative of your report, but the report should be in the style of a self containing paper, not a numbered list to the questions here.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b944b78b",
   "metadata": {},
   "source": [
    "(1) Introduction\n",
    "----------------\n",
    "a) Read the paper by Allison (2012) and discuss with your peers. What did the author do?  What was the motivation for it? Which terms or concept are new to you?  Look them up.\n",
    "\n",
    "b) What is the Gaia mission? What will/ has Gaia measured?  What Gaia measurements are relevant for this paper?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9eadb45a",
   "metadata": {},
   "source": [
    "(2) Initial conditions\n",
    "-----------------------\n",
    "\n",
    "a) Lookup the references given in the paper of the initial conditions for the spatial distribution of stars and their mass distribution. Keep them for reference.\n",
    "\n",
    "Luckily the initial conditions are already included with AMUSE!\n",
    "\n",
    "b) Import the fractal model cluster, and read the help to get information on the arguments of this function:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80371773",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from amuse.ic.fractalcluster import new_fractal_cluster_model\n",
    "help(new_fractal_cluster_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ae6ea97",
   "metadata": {},
   "source": [
    "c) Generate a N=100 model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d75fb01",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "parts=new_fractal_cluster_model(100, do_scale=True)\n",
    "print(parts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a75068b3",
   "metadata": {},
   "source": [
    "d) Do you remember from the tutorial what the units are? Confirm that the total mass and potential and kinetic energies match the expected values. (The potential energy calculation is given below - what happens if you do not provide G? can you understand the generated error?)\n",
    "\n",
    "e) explore the effect of the fractal dimension and virial ratio parameters on the model by plotting the positions and \n",
    "velocities of various models. Do you understand the physical meaning of these?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d6c838f",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from amuse.units import units,constants, nbody_system\n",
    "\n",
    "parts.potential_energy(G=nbody_system.G)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f780a3a5",
   "metadata": {},
   "source": [
    "To finish the initial conditions we need to convert the model to physical units, and give the stars a realistic distribution of masses.\n",
    "\n",
    "f) The Kroupa mass distribution used in the paper can be generated by a a piece-wise power law distribution. Generate a set of masses described by this, and plot the resulting mass distribution.Adapt the parameters to match the paper!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc02d872",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from amuse.ic.brokenimf import new_broken_power_law_mass_distribution\n",
    "masses=new_broken_power_law_mass_distribution(1000, mass_boundaries=[0.1,100] | units.MSun, alphas=[-2.35])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d248b99",
   "metadata": {},
   "source": [
    "g) The way to convert a scaleless model to SI units is by providing a converter. An N-body converter takes two values, fixing the scale of the model (together with the requirement of G=1). (can you understand the units of the converted values below?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8aaf878a",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "radius=1 | units.parsec\n",
    "\n",
    "convert=nbody_system.nbody_to_si( masses.sum(), radius)\n",
    "\n",
    "print(convert.to_si( 1 | nbody_system.length ))\n",
    "print(convert.to_nbody( 1 | units.km/units.s ))\n",
    "print(convert.to_nbody( constants.G ))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ef9e6ae",
   "metadata": {},
   "source": [
    "h) at this point you should be able to generate intitial conditions matching those of the paper! Do so."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d898b1bf",
   "metadata": {},
   "source": [
    "(3) Evolving the models\n",
    "-----------------------\n",
    "\n",
    "a) Like the paper mentions, the experiments where done using purely gravitational dynamics. A number of codes for this \n",
    "are available in AMUSE. Choose an appropiate one from the tabel given in the AMUSE code paper. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c5c6bcd",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from amuse.community.huayno.interface import Huayno\n",
    "\n",
    "h=Huayno(convert_nbody=convert)\n",
    "print(h.parameters)\n",
    "help(h.parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c72d066",
   "metadata": {},
   "source": [
    "b) write an evolve loop to evolve the model to some specified end_time with some dt, measuring and reporting the energy error \n",
    "\n",
    "c) conduct the simulations as reported in the paper. It may be that a simulation takes too long or the energy error is not acceptable ((E-E0)/E0 should be smaller than 0.001). You can try with smaller number of particles, introducing smoothing (epsilon_squared parameter) or change the timestep parameter (usually timestep_parameter), or vary the random seed of the initial distribution (can you figure out the rationale behind these suggestions?)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8715c766",
   "metadata": {},
   "source": [
    "(4) Analysis\n",
    "------------\n",
    "\n",
    "a) reproduce the spatial distribution plot of Allison 2012. Note similarities and differences\n",
    "\n",
    "b) analyse the kinematics of the escaping stars, discussing Allision 2012.\n",
    "\n",
    "Explore a bit of parameter space if you feel like it!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dc08553",
   "metadata": {},
   "source": [
    "(5) Wrapping it up\n",
    "------------------\n",
    "\n",
    "a) This exercise is completed by writing short report, see instructions below.  You can distill from the above notebook (maybe copying from the above to a new notebook). Include 1 or 2 figures for the results (matching those from the paper) and maybe 1 or 2 more showing the initial conditions.\n",
    "\n",
    "b) If you didn't succeed, try to describe why not and what the problems where you ran into. \n",
    "\n",
    "c) the orginal paper came out well before the GAIA data releases. Try to find out if anyone has followed up the work with the currently available GAIA data.\n",
    "\n",
    "d) Bonus points if you extend the model with any of the additional physics mentioned in the paper (stellar evolution, binaries natal gas, but in particular the galactic potential)! But first make sure you full fill the basics. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a58f078",
   "metadata": {},
   "source": [
    "Format of your report\n",
    "---------------------\n",
    "\n",
    "Write your report as a \"letter\" (= a brief paper of 4-5 pages max in 2 column format when using the  ``emulateapj`` package. This is a latex class that will give it the appearance of a real ApJ letter. You can download it online. I suggest you use the following latex packages\n",
    "```\n",
    "    \\\\documentclass[apjl]{emulateapj}  % paper will get the appearance of a real ApJ Letter\n",
    "    \\\\usepackage{graphicx} % for including plots\n",
    "    \\\\usepackage{natbib} % for a proper bibiography\n",
    "    \\\\bibliographystyle{apj} % to set the style of the bibliography\n",
    "    \\\\usepackage{hyperref} % will make all references in your paper act as hyperlinks in the pdf\n",
    "```\n",
    "Make sure you have the following sections and do not exceed 5 pages.\n",
    "    1. Introduction\n",
    "    2. Method\n",
    "    3. Results\n",
    "    4. Conclusion & Discussion.\n",
    "\n",
    "Some hints for writing\n",
    "----------------------\n",
    "\n",
    "##### Introduction:\n",
    "A good introduction briefly explains\n",
    "\n",
    "* the astrophysical context / the bigger question your work adresses / why your work is timely (hint: maybe because it relevant for the data that Gaia will soon release)\n",
    "* the state of the art and gives credit to most relevant earlier papers (Hint: it would be nice to state that you follow the approach taken by Allison 2012)\n",
    "* formulate the specific scientific question you are adressing in this work.\n",
    "* is concise: aim for 0.5-0.75 page in the 2 col format\n",
    "\n",
    "##### Method\n",
    "\n",
    "Here you explain which code(s) you are using and you give a very brief summary of what they do providing references. Then state the assumptions and initial conditions you adopt. Typical length is ~1-1.5 page in a letter.\n",
    "\n",
    "##### Results\n",
    "\n",
    "Here you describe your main results by systematically dicussing your findings. Hint: It often helps to first decide what the key plots are that you want to show.  Often the result section is often a systematic explanationof what the plots illustate.  Typical length is 2 pages in a letter.\n",
    "\n",
    "##### Discussion & conclusion \n",
    "Summarize your key message. Here you can disuss the implications of your findings, the limitations of your method and possibly state future improvements. Typical length is 0.5 pages in a letter.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68305997",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "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": 5
}