10 - Generic_units.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=\"amuse-framework 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": 1,
   "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": [
    "## Generic and Nbody unit systems\n",
    "\n",
    "In this tutorial the generic and nbody systems are explained in more detail"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the tutorial on units and quantities you have familiarized yourself with the use of units and the construction of quantities which are a combination of numerical values and an associated unit. The use of units in AMUSE serves two purposes: 1) it gives physical meaning to the results of simulation codes and the outcome of calculations and 2) they serves mechanism to check the consistency of the calculations. The physical meaning allows us two codes in AMUSE to understand each others output (e.g. conversion from m to parsec), the consistency check makes sure that the input are meaningful (ie no masses are used as positions).\n",
    "\n",
    "In many community codes and for a lot of algorithms it makes sense to use units for which the base units are not (fully) specified. This is often the case when e.g. the equations solved are scale free or for initial conditions where similar models can be scaled to different sizes. Although the quantities involved in this case do not have a specific unit base, they still have a dimension (mass, length etc.). So in this case, the variables of a code, e.g. the position, do not have a particular scale, they still have the dimension length. In AMUSE we can use *generic units* in this case. In other words, you can specify if a value has a *mass*, *length* or *time* dimension, or any combination thereof, such as *length* per *time*. \n",
    "\n",
    "We will make this clearer with a few examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.units import units, constants, nbody_system, generic_unit_system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "AMUSE includes two generic unit systems, the **generic_unit_system** is the most general, the **nbody_system** is a special case and defines the gravitational constant to have the numerical value `G=1` (as we will see G in this case still has dimensions!). For gravity calculations the **nbody_system** module is recommended as this follows the general practice in most n-body codes. Note that in the literature N-body units often have the additional property that he system being studied is scaled such that total mass `M=1` and energy `E=-0.25` (there are also called Hénon units).\n",
    "\n",
    "The generic units are defined in the **generic_system** and **nbody_system** modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(10.0 | nbody_system.length)\n",
    "print(10.0 | generic_unit_system.length)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quantities with generic units work exactly the same as quantities with normal (**S.I.**) units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_mass = 1.0 | generic_unit_system.length\n",
    "mean_speed = 0.1 | generic_unit_system.length / generic_unit_system.time\n",
    "print(mean_speed * cluster_mass)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Generic quantities are very useful and can be applied almost everywhere in AMUSE. \n",
    "\n",
    "To convert quantities to and from a specific system of units you'll need a converter. For nbody units you can create a converter to S.I. units like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "converter = nbody_system.nbody_to_si(1 | units.MSun, 1 | units.kms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An ``nbody_system`` converter is specified by passing to it two independent quantities, which are to be scaled to 1. Together with the requirement that `G=1` this fixes the unit system. These quantities can be simple (like 1 solar mass) or combined (like  10 km/s). With the scaling the resulting object converts to and from the nbody units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Mass of the sun, scaled:\", converter.to_nbody(1 | units.MSun))\n",
    "print(\"10 nbody masses, in S.I.:\", converter.to_si(10 | nbody_system.mass))\n",
    "print(\"1 nbody time, in S.I:\", converter.to_si(1 | nbody_system.time).in_(units.yr))\n",
    "print(\"10 km/s, in nbody:\", converter.to_nbody(10.0 | units.km / units.s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets check the scaling explicitly: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"For {converter}:\")\n",
    "print(f\" G = {converter.to_nbody(constants.G)}\")\n",
    "print(f\" 1 MSun = {converter.to_nbody(1 | units.MSun)}\")\n",
    "print(f\" 1 km/s = {converter.to_nbody(1| units.kms)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the nbody units of these quantities all have the numerical value `1` (however they are not dimensionless! scale-free does not mean dimensionless!!)\n",
    "\n",
    "It is important to understand what the converter does: it provides a particular scaling from a generic unit system to (most often) S.I. . In AMUSE we use it often in two roles: passed as an argument to a code class constructor it fixes the scaling of the input and output quantities of the code (be it variables or parameters). Passed as an argument to an initial condition generator it fixes the scaling of that initial condition generator. So this is very useful - but it is not strictly required! (we will give examples further down)   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the generic unit converter, you can specify up to 7 quantities (as there are 7 base properties). Any combination of quantities is possible as long as it results in a orthogonal set of converters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "converter = generic_unit_converter.ConvertBetweenGenericAndSiUnits(1 | units.MSun, 1 | units.AU, constants.G)\n",
    "print(\"Mass of the sun, scaled:\", converter.to_nbody(1 | units.MSun))\n",
    "print(\"10 generic masses, in S.I.:\", converter.to_si(10 | nbody_system.mass))\n",
    "print(\"1 generic time, in S.I:\", converter.to_si(1 | nbody_system.time).in_(units.yr))\n",
    "print(\"10 km/s, in generic:\", converter.to_nbody(10.0 | units.km / units.s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "converter = generic_unit_converter.ConvertBetweenGenericAndSiUnits(1 | units.MSun, 1 | units.AU, 1 | units.yr)\n",
    "print(\"Mass of the sun, scaled:\", converter.to_nbody(1 | units.MSun))\n",
    "print(\"10 generic masses, in S.I.:\", converter.to_si(10 | nbody_system.mass))\n",
    "print(\"1 generic time, in S.I:\", converter.to_si(1 | nbody_system.time).in_(units.yr))\n",
    "print(\"10 km/s, in generic:\", converter.to_nbody(10.0 | units.km / units.s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specifying a length twice or specifying a speed and a length and a time will result in an error. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "generic_unit_converter.ConvertBetweenGenericAndSiUnits(\n",
    "    1 | units.MSun, \n",
    "    1 | units.AU, \n",
    "    1 | units.m,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, the following defines a converter for Planck units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "natural_units_convert = generic_unit_converter.ConvertBetweenGenericAndSiUnits(\n",
    "    constants.c,\n",
    "    constants.G,\n",
    "    constants.hbar,\n",
    "    1/(4*numpy.pi*constants.eps0),\n",
    "    constants.kB,\n",
    ")\n",
    "\n",
    "M = 1 | generic_unit_system.mass\n",
    "T = 1 | generic_unit_system.time\n",
    "L = 1 | generic_unit_system.length\n",
    "Q = 1 | generic_unit_system.charge\n",
    "THETA = 1 | generic_unit_system.temperature\n",
    "\n",
    "print(natural_units_convert.to_si(M).in_(units.kg))\n",
    "print(natural_units_convert.to_si(T).in_(units.s))\n",
    "print(natural_units_convert.to_si(L).in_(units.m))\n",
    "print(natural_units_convert.to_si(Q).in_(units.C))\n",
    "print(natural_units_convert.to_si(THETA).in_(units.K))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course unit commensurability is still enforced:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print((10.0 | nbody_system.length) + (10.0 | generic_unit_system.time))"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 1\n",
    "\n",
    "Examine and try out the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = new_plummer_model(100)\n",
    "h = Huayno()\n",
    "print(h.parameters)\n",
    "h.particles.add_particles(p)\n",
    "h.evolve_model(1 | nbody_system.time)\n",
    "print(h.particles)\n",
    "h.stop()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "convert = nbody_system.nbody_to_si(100. | units.MSun, 1. | units.parsec)\n",
    "p = new_plummer_model(100, convert_nbody=convert)\n",
    "h = Huayno(convert)\n",
    "print(f\"smoothing parameter: {h.parameters.epsilon_squared}\")\n",
    "h.particles.add_particles(p)\n",
    "h.evolve_model(1| units.Myr)\n",
    "print(h.particles)\n",
    "h.stop()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "convert = nbody_system.nbody_to_si(100. | units.MSun, 1. | units.parsec)\n",
    "p = new_plummer_model(100)\n",
    "h = Huayno(convert)\n",
    "print(f\"smoothing parameter: {h.parameters.epsilon_squared}\")\n",
    "h.particles.add_particles(p)\n",
    "h.evolve_model(1| nbody_system.time)\n",
    "print(h.particles)\n",
    "print(f\" reached {h.model_time.in_(units.Myr)}\")\n",
    "h.stop()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 2\n",
    "\n",
    "Look up and listen to the Huayno preserved on the Voyager 2 Golden record."
   ]
  },
  {
   "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
}