05 - Running_stellar_evolution.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "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\"\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": {
    "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 stellar evolution code\n",
    "====================\n",
    "\n",
    "Here we will generate initial conditions for a stellar evolution code. We will adopt a parametrized stellar evolution code.\n",
    "\n",
    "For reference you can read chapter 3 of Portegies Zwart & McMillan 2018 (2018araa.book.....P). \n",
    "\n",
    "With this tutorial you will learn\n",
    "- how to generate inital conditions (mass-function)\n",
    "- how to run a parametrized stellar evolution code\n",
    "- channels and data transfer\n",
    "- plotting results"
   ]
  },
  {
   "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\n",
    "from amuse.datamodel import Particles\n",
    "from amuse.community.seba.interface import Seba\n",
    "from amuse.community.sse.interface import Sse\n",
    "from amuse.lab import new_kroupa_mass_distribution\n",
    "from amuse.lab import new_salpeter_mass_distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This sets up and evolves a 1 solar mass star:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p=Particles(mass=[1] | units.MSun)\n",
    "s=Sse()\n",
    "p=s.particles.add_particles(p)\n",
    "s.evolve_model(1| units.Myr)\n",
    "print(p.age, p.luminosity, p.temperature,p.time_step)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 1:\n",
    "\n",
    "write a code snippet that plots a Hertzsprung-Russel diagram (log-log plot of luminosity vs temperature) for the evolutionary track of a star of given mass. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stellar populations \n",
    "\n",
    "Ofcourse, we may be interested in the behaviour of stellar populations. Lets work through a more interesting example.\n",
    "\n",
    "We start by setting up a simulation by specifying the initial mass function. What would you say is the biggest difference between these IMFs?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_stars = 1024\n",
    "mmin = 0.1 | units.MSun\n",
    "mmax = 100 | units.MSun\n",
    "\n",
    "mkroupa = new_kroupa_mass_distribution(\n",
    "    n_stars,\n",
    "    mass_min=mmin, \n",
    "    mass_max=mmax)\n",
    "k_stars = Particles(mass=mkroupa)\n",
    "\n",
    "msalpeter = new_salpeter_mass_distribution(\n",
    "    n_stars,\n",
    "    mass_min=mmin,\n",
    "    mass_max=mmax)\n",
    "s_stars = Particles(mass=msalpeter)\n",
    "print(\"mean-mass for Kroupa=\", numpy.mean(k_stars.mass))\n",
    "print(\"mean-mass for Salpeter=\", numpy.mean(s_stars.mass))\n",
    "\n",
    "pyplot.hist(numpy.log10(k_stars.mass.value_in(units.MSun)))\n",
    "pyplot.title('Kroupa IMF')\n",
    "pyplot.yscale('log')\n",
    "pyplot.xlim(-1, 2)\n",
    "pyplot.figure()\n",
    "pyplot.hist(numpy.log10(s_stars.mass.value_in(units.MSun)))\n",
    "pyplot.title('Salpeter IMF')\n",
    "pyplot.yscale('log')\n",
    "pyplot.xlim(-1, 2)\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the mass functions are determined, we can evolve the stars in time. For this we start two stellar evolution codes, both parametrized, and evolve the stars in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "se = Sse()\n",
    "kstellar = se.particles.add_particles(k_stars)\n",
    "sstellar = se.particles.add_particles(s_stars)\n",
    "\n",
    "kchannels = {\"to_stars\": kstellar.new_channel_to(k_stars), \n",
    "                \"to_stellar\": k_stars.new_channel_to(kstellar)}\n",
    "schannels = {\"to_stars\": sstellar.new_channel_to(s_stars), \n",
    "                \"to_stellar\": s_stars.new_channel_to(sstellar)}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the stellar evolution code is declared. Now evolve in time and calculate the mean mass as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = 10**numpy.arange(0.0, 4.0, 0.1) | units.Myr\n",
    "mmean = []\n",
    "for time in times:\n",
    "    print(f\"evolving to {time}\")\n",
    "    se.evolve_model(time)\n",
    "    kchannels[\"to_stars\"].copy()\n",
    "    kchannels[\"to_stars\"].copy_attributes([\"temperature\"]) # needs explicit copy with SSE because its derived\n",
    "    schannels[\"to_stars\"].copy()\n",
    "    schannels[\"to_stars\"].copy_attributes([\"temperature\"])\n",
    "    mmean.append(numpy.mean(k_stars.mass)/numpy.mean(s_stars.mass))\n",
    "se.stop()\n",
    "print(\"mean-mass for Kroupa=\", numpy.mean(k_stars.mass))\n",
    "print(\"mean-mass for Salpeter=\", numpy.mean(s_stars.mass))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.plot import plot\n",
    "plot(times, mmean)\n",
    "pyplot.ylabel(\"relative mean mass\")\n",
    "pyplot.semilogx()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.plot import scatter\n",
    "scatter(s_stars.temperature, s_stars.luminosity, c='r')\n",
    "scatter(k_stars.temperature, k_stars.luminosity, c='b', s=3)\n",
    "pyplot.xlim(2.e+4, 2000)\n",
    "pyplot.ylim(1.e-5, 1000)\n",
    "pyplot.loglog()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have calculated the evolution of a stellar population with solar metalicity.\n",
    "\n",
    "Assignments and questions:\n",
    "---------------\n",
    "\n",
    "### Assignment 2:\n",
    "Run the notebook and describe the evolution of the relative mean stellar mass. Why does it drop with time?\n",
    "\n",
    "### Assignment 3:\n",
    "Compare the mean mass of two stellar populations with different metalicity, 0.02 (solar) and 0.002 (one tenth times solar) but the same stellar masses.\n",
    "\n",
    "### Question 1:\n",
    "Is there a difference in the number of compact objects (white dwarfs, neutron stars, black holes; these types are defined [here](https://github.com/amusecode/amuse/blob/main/src/amuse/units/units.py), and e.g. 1.|stellar_type corresponds to a main sequene star) in the populations with different initial mass functions, and with different metallicities, and is this difference statistically significant?\n",
    "\n",
    "Is a similar difference also noticeable in the number of main sequence stars?\n",
    "\n",
    "### Assignment 4:\n",
    "Calculate, instead of the mean mass, the total luminosity and luminosity-weighted temperature ($T_{LW}=(\\Sigma_iL_iT_i)/(\\Sigma_iL_i)$, the effective color of the integrated stellar population) of the entire population of stars as a function of time.\n",
    "\n",
    "Make a plot of these parameters for both mass functions, and both metallicities. \n",
    "\n",
    "### Question 4:\n",
    "Can one, based on observational differences in luminosity and integrated temperature, make a distinction between the two mass functions?"
   ]
  },
  {
   "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
}