high_order_bridge.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "AMUSE tutorial on high-order bridge\n",
    "====================\n",
    "\n",
    "Hierarchical coupling strategies are fundamental parts of AMUSE.\n",
    "It enables us to combine the output of a wide variety of solvers into a homogeneous solution.\n",
    "In this example we will be nesting multiple bridges, to show the power of bridge.\n",
    "\n",
    "We realize this by integrating a star that is orbited by a planet, which is orbited by a moon.\n",
    "In principle, of course, one can put them all in a single N-body integratory, but here we use bridge.\n",
    "\n",
    "Let's start by importing the relevant packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "from amuse.units import (units, constants)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a convenient function for obtaining the orbital period of a binary system, and initilize a system with star, planet and moon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def orbital_period(Mtot, a):\n",
    "    return (((4 * numpy.pi**2) * a**3)/(constants.G * Mtot)).sqrt()\n",
    "\n",
    "def Hill_radius(Mstar, Mplanet, a_planet, e_planet):\n",
    "    return a_planet*(1.0-e_planet)*(Mplanet/(3*Mstar))**(1./3.)\n",
    "    \n",
    "from amuse.ext.orbital_elements import new_binary_from_orbital_elements\n",
    "def get_star_planet_and_moon():\n",
    "    Mstar = 1.0|units.MSun\n",
    "    Mplanet = 1.0|units.MJupiter\n",
    "    a_planet = 5.2 | units.au\n",
    "    e_planet = 0.0\n",
    "    bodies = new_binary_from_orbital_elements(Mstar, Mplanet, a_planet, e_planet,\n",
    "                                             G=constants.G)\n",
    "    bodies[0].name = \"star\"\n",
    "    bodies[1].name = \"planet\"\n",
    "    planet = bodies[bodies.name==\"planet\"]\n",
    "    RH_planet = Hill_radius(Mstar, Mplanet, a_planet, e_planet)\n",
    "\n",
    "    a_moon = 0.1*RH_planet\n",
    "    e_moon = 0.0\n",
    "    Mmoon = 0.01*Mplanet\n",
    "    Pmoon = orbital_period(Mplanet+Mmoon, a_moon)\n",
    "    moon = new_binary_from_orbital_elements(planet.mass, \n",
    "                                            Mmoon, a_moon, e_moon,\n",
    "                                            G=constants.G)\n",
    "    moon.position -= moon[0].position\n",
    "    moon.velocity -= moon[0].velocity\n",
    "    moon = moon[1].as_set()\n",
    "    moon.position += planet.position\n",
    "    moon.velocity += planet.velocity\n",
    "    moon.semimajor_axis = a_moon\n",
    "    moon.name = \"moon\"\n",
    "    bodies.add_particle(moon)\n",
    "    bodies.move_to_center()\n",
    "    return bodies\n",
    "starplanmoon = get_star_planet_and_moon()\n",
    "print(starplanmoon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the orbits of the three particles, star, planet and moon. to make things a bit more interesting, we make a debris disk around the moon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.ext.protodisk import ProtoPlanetaryDisk\n",
    "\n",
    "from amuse.lab import nbody_system\n",
    "\n",
    "def make_disk_around_celestial_body(moon, Mplanet):\n",
    "\n",
    "    a_moon = moon.semimajor_axis[0]\n",
    "    e_moon = 0.0\n",
    "    Mmoon = moon.mass.sum()\n",
    "    RH_moon = Hill_radius(Mplanet, Mmoon, a_moon, e_moon)\n",
    "    converter = nbody_system.nbody_to_si(Mmoon, RH_moon)\n",
    "    Ndisk = 100\n",
    "    Rin = 0.03\n",
    "    Rout = 1.0\n",
    "    Pinner = orbital_period(Mmoon, Rin * a_moon)\n",
    "    Mdisk = 0.01\n",
    "\n",
    "    disk = ProtoPlanetaryDisk(Ndisk,\n",
    "                              convert_nbody=converter,\n",
    "                              Rmin=Rin,\n",
    "                              Rmax=Rout,\n",
    "                              q_out=10.0,\n",
    "                              discfraction=Mdisk).result\n",
    "    disk.move_to_center()\n",
    "    disk.position += moon.position\n",
    "    disk.velocity += moon.velocity\n",
    "\n",
    "    disk.mass = Mdisk*Mmoon/float(Ndisk)\n",
    "    rho = 3.0| (units.g/units.cm**3)\n",
    "    disk.radius = (disk.mass/(4*rho))**(1./3.)\n",
    "    disk.name = \"disk\"\n",
    "    return disk, Pinner, converter\n",
    "\n",
    "planet = starplanmoon[starplanmoon.name==\"planet\"]\n",
    "moon = starplanmoon[starplanmoon.name==\"moon\"]\n",
    "disk, Pinner, converter = make_disk_around_celestial_body(moon, \n",
    "                                               planet.mass.sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just define some handy routines for plotting and printing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot\n",
    "def plot_system(bodies):\n",
    "    star = bodies[bodies.name==\"star\"]\n",
    "    planets = bodies[bodies.name==\"planet\"]\n",
    "    moon = bodies[bodies.name==\"moon\"]\n",
    "    disk = bodies[bodies.name==\"disk\"]\n",
    "    pyplot.scatter(star.x.value_in(units.au), star.y.value_in(units.au), c='y', s=100)\n",
    "    pyplot.scatter(planet.x.value_in(units.au), planet.y.value_in(units.au), c='b', s=30)\n",
    "    pyplot.scatter(disk.x.value_in(units.au), disk.y.value_in(units.au), c='k', s=1)\n",
    "    pyplot.scatter(moon.x.value_in(units.au), moon.y.value_in(units.au), c='r', s=10)\n",
    "    pyplot.axis(\"equal\")\n",
    "    pyplot.show()\n",
    "    \n",
    "def print_system(bodies):\n",
    "    star = bodies[bodies.name==\"star\"]\n",
    "    planets = bodies[bodies.name==\"planet\"]\n",
    "    moon = bodies[bodies.name==\"moon\"]\n",
    "    disk = bodies[bodies.name==\"disk\"]\n",
    "    print(\"moon:\", moon.position.in_(units.au))\n",
    "    print(\"disk:\", disk.center_of_mass().in_(units.au))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Declare the various gravity solvers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.community.huayno.interface import Huayno\n",
    "bodies = starplanmoon\n",
    "bodies.add_particles(disk)\n",
    "\n",
    "star = bodies[bodies.name==\"star\"]\n",
    "planets = bodies[bodies.name==\"planet\"]\n",
    "moon = bodies[bodies.name==\"moon\"]\n",
    "disk = bodies[bodies.name==\"disk\"]\n",
    "    \n",
    "gravityA = Huayno(converter)\n",
    "gravityA.particles.add_particles(star)\n",
    "gravityA.particles.add_particles(planets)\n",
    "channel = {\"from stars\": bodies.new_channel_to(gravityA.particles),\n",
    "            \"to_stars\": gravityA.particles.new_channel_to(bodies)}\n",
    "\n",
    "gravityB = Huayno(converter)\n",
    "gravityB.particles.add_particles(moon)\n",
    "gravityB.particles.add_particles(disk)\n",
    "channel.update({\"from_disk\": bodies.new_channel_to(gravityB.particles)})\n",
    "channel.update({\"to_disk\": gravityB.particles.new_channel_to(bodies)})\n",
    "\n",
    "plot_system(bodies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Declare the bi-directional bridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.couple import bridge\n",
    "from amuse.ext.composition_methods import *\n",
    "gravity = bridge.Bridge(use_threading=False)\n",
    "gravity.add_system(gravityA, (gravityB,))\n",
    "gravity.add_system(gravityB, (gravityA,))\n",
    "gravity.timestep = 0.1*Pinner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the code in the event loop. Store some data and print some diagnostics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.ext.composition_methods import *\n",
    "from amuse.ext.orbital_elements import orbital_elements_from_binary\n",
    "\n",
    "def gravity_gravity_bridge(gravityA, gravityB, gravity, bodies,\n",
    "                         t_end):\n",
    "\n",
    "    gravity_initial_total_energy = gravityA.get_total_energy() + gravityB.get_total_energy()\n",
    "    model_time = 0 | units.Myr\n",
    "    dt = t_end/20.\n",
    "    \n",
    "    while model_time < t_end:\n",
    "\n",
    "        model_time += dt\n",
    "\n",
    "        orbit_planet = orbital_elements_from_binary(bodies[:2], G=constants.G)\n",
    "        orbit_moon = orbital_elements_from_binary(bodies[1:3], G=constants.G)\n",
    "        print(\"Planet:\", \"a=\", orbit_planet[2].in_(units.AU), \"e=\", orbit_planet[3])\n",
    "        print(\"Moon:\", \"a=\", orbit_moon[2].in_(units.AU), \"e=\", orbit_moon[3])\n",
    "        \n",
    "        dE_gravity = gravity_initial_total_energy/(gravityA.get_total_energy()+gravityB.get_total_energy())\n",
    "        print(\"Time:\", model_time.in_(units.day), \\\n",
    "              \"dE=\", dE_gravity)#, dE_hydro\n",
    "\n",
    "        gravity.evolve_model(model_time)\n",
    "        channel[\"to_stars\"].copy()\n",
    "        channel[\"to_disk\"].copy() \n",
    "        print_system(bodies)\n",
    "        print(\"com=\", gravityB.particles.center_of_mass().in_(units.au))\n",
    "        plot_system(bodies)\n",
    "    \n",
    "t_end = 100.0 | units.day\n",
    "gravity_gravity_bridge(gravityA, gravityB, gravity, \n",
    "                       bodies, t_end)\n",
    "gravity.stop()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have created a debris around a moon in orbit around a planet that orbits a star.\n",
    "In principle, such systems have not yet been observed, but we can simulate them anyway.\n",
    "The code you have been running can be adapted and improved.\n",
    "\n",
    "\n",
    "Assignments and questions:\n",
    "---------------\n",
    "\n",
    "Huayno is usually 2nd order. But we can at least increase the order of the coupling between the codes in the bridge.\n",
    "This is done using the methods keyword in the bridge operator.\n",
    "Integrate the system using the SPLIT_4TH_S_M4 or SPLIT_10TH_SS_M35 method (or one of the other methods). These two are 4th and 10th order integration scheme within the bridge.\n",
    "One advantage of a higher-order scheme is the possibility to increase the time step, and still get a pretty good results. The quality of the result is generally indicated with the degree by which energy is conserved, because we consider energy conservation an essential quality of the reliable numerical solution.\n",
    "Increasing the order of the integrator then allows us to increase the time step while keeping the error in the energy the same. The integration scheme of the standard bridge is second order. In the above experiment we adopted a time step of 0.001 of the orbital period of the moon. By adopting a higher order bridge coupling we can increase this order.\n",
    "\n",
    "### Assignment 1:\n",
    "The plotting routine is somewhat simple. Rewrite the routine\n",
    "plot_system in such a way that you can study the disk in more detail.\n",
    "This can be achieved by centering the picture around the planet or on the moon.\n",
    "\n",
    "Another interesting diagnostic would be the evolution of the moon's orbital parameters.\n",
    "Make a plot of the moon's semi-major axis and eccentricity in time.\n",
    "\n",
    "### Assignment 2:\n",
    "\n",
    "The disk seems to be rather unstable, and the majority of material evaporates within a few orbits.\n",
    "Adjust the disk in such a way, that it becomes more stable.\n",
    "\n",
    "### Question 1:\n",
    "\n",
    "How small do you have to make the disk before it is stable for one orbit?\n",
    "Does it help to make the time step smaller to acquire a more stable disk?\n",
    "\n",
    "\n",
    "### Question 2:\n",
    "\n",
    "The code is not fast.  What part of the code is slow and why?\n",
    "One way to make the code faster is by increasing the order of the bridge integrator while increasing the time step.\n",
    "Adopt a 4th order bridge integrator with a time step of 0.01 of the inner orbital period. What effect does this has on the performance and energy conservation when running the code?\n",
    "\n",
    "Do the same but then with a 10th order integration scheme in the bridge and a bridge time step of 0.1 of the Moon's orbit. What affect that this has on performance and accuracy?\n",
    "\n",
    "### Assignment 3:\n",
    "\n",
    "You can improve the code's performance in three ways.\n",
    "\n",
    " * One is by running Huayno in parallel. This is done via openMP, not using MPI as most packages in AMUSE. Time the codes' performance, and then run with openMP. \n",
    "\n",
    " * Another way to improve Huayno's performance is by calculating the orbital evolution of the test particles using a Kepler solver rather than integrating the orbit. The bridge will take care of the interaction between disk particles and the planet. The name of the solver is CC_KEPLER, indicating that it adopts a connected-components strategy while solving two-particle interactions using Kepler orbits. \n",
    " \n",
    "You can change the type of integrator in Huayno using:\n",
    "gravityB.parameters.inttype_parameter = Huayno.inttypes.CC_KEPLER\n",
    "\n",
    " * Eventually, speed-up can be achieved by making the particles in the disk zero-mass. In Huayno such particles are automatically recognized as test particles and not integrated with the N^2 algorithm.\n",
    "\n",
    "### Question 3:\n",
    "\n",
    "When adopting a different way of integrating your system, either by making the debris particles zero-mass, or by introducing the CC_KEPLER solver will affect your results. This is noticeable in the morphology of the disk and the number of escaping particles. How can you still validate your results if the outcome depends on the solver?\n",
    "\n",
    "How much faster do you expect the code is when running with zero-mass disk particles rather than with finite mass?\n",
    "Why does using openMP not affect your results?\n",
    "\n",
    "### Assignment 4:\n",
    "\n",
    "The planetary system we are studying may have a secondary moon. Add a second moon to the planet in a circular orbit at a distance of 0.3 Hill radii around the planet. \n",
    "\n",
    "### Question 4:\n",
    "\n",
    "Where do you put this secondary moon, in the gravityA, in the gravityB or do you bridge with yet another gravity solver?\n",
    "\n",
    "### Question 5:\n",
    "\n",
    "What effect does this moon have on the debris disk around the first moon?\n",
    "\n",
    "### Assignment 5:\n",
    "\n",
    "Now the star may actually have an equal-mass companion star. You want to keep the orbit of your planet stable. If this secondary star has a circular orbit in the plane of the planet, at what distance should it orbit to prevent it from perturbing the planet?\n",
    "\n",
    "### Question 6: \n",
    "\n",
    "How is the disk around the moon in orbit around the planet affected by the presence of this secondary star?\n",
    "\n",
    "### Assignment 6:\n",
    "\n",
    "Debris may fall onto one of the two moons or onto the planet.\n",
    "Add a stopping condition for collision detection for the star, planet and moons.\n",
    "Assume all objects to have the same density of 1g/cc.\n",
    "\n",
    "### Question 7:\n",
    "\n",
    "Note that it turns out to be tricky to introduce a stopping conditions between two particle sets that live in other codes. How can you solve this problem?\n",
    "\n",
    "Which of the objects, star, planet or which moon, is most frequenly hit by the debris? \n",
    "Can you explain this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}