bridge_gravity_with_hydro.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 combined 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "from amuse.units import (units, constants)\n",
    "#from amuse.lab import Particles"
   ]
  },
  {
   "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",
    "from amuse.ext.orbital_elements import new_binary_from_orbital_elements\n",
    "from amuse.ext.orbital_elements import orbital_elements_from_binary\n",
    "\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 = a_planet * (1.0 - e_planet) * (Mplanet / (3 * Mstar))**(1./3.)\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",
    "bodies = get_star_planet_and_moon()\n",
    "print(bodies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the orbits of the three particles, star, planet and moon. We now want to make a 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",
    "\n",
    "def make_disk_around_celestial_body(moon, Mplanet):\n",
    "\n",
    "    R = 1 | units.au\n",
    "    a_moon = moon.semimajor_axis\n",
    "    e_moon = 0.0\n",
    "    Mmoon = moon.mass.sum()\n",
    "    RH_moon = a_moon * (1.0 - e_moon) * (Mmoon / (3 * Mplanet))**(1./3.)\n",
    "    converter = nbody_system.nbody_to_si(Mmoon, R)\n",
    "    Ndisk = 1000\n",
    "    Rin = 0.03 * RH_moon\n",
    "    Rout = 0.3 * RH_moon\n",
    "    Pinner = orbital_period(Mmoon, Rin)\n",
    "    Mdisk = 0.01 * Mmoon\n",
    "\n",
    "    disk = ProtoPlanetaryDisk(Ndisk,\n",
    "                              convert_nbody=converter,\n",
    "                              Rmin=Rin/R,\n",
    "                              Rmax=Rout/R,\n",
    "                              q_out=10.0,\n",
    "                              discfraction=Mdisk/Mmoon).result\n",
    "    disk.name = \"disk \"\n",
    "    disk.move_to_center()\n",
    "    disk.position += moon.position\n",
    "    disk.velocity += moon.velocity\n",
    "\n",
    "    masses = Mdisk / float(Ndisk)\n",
    "    disk.mass = masses\n",
    "    rho = 3.0 | (units.g / units.cm**3)\n",
    "    disk.radius = (disk.mass / (4 * rho))**(1./3.)\n",
    "    return disk, Pinner, converter\n",
    "\n",
    "planet = bodies[bodies.name==\"planet\"]\n",
    "moon = bodies[bodies.name==\"moon\"]\n",
    "disk, Pinner, converter = make_disk_around_celestial_body(moon, \n",
    "                                               planet.mass.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.plot import scatter\n",
    "def plot(moon, disk):\n",
    "    scatter(disk.x, disk.y)\n",
    "    scatter(moon.x, moon.y, c='r')\n",
    "plot(moon, disk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.community.huayno.interface import Huayno\n",
    "from amuse.community.fi.interface import Fi\n",
    "gravity = Huayno(converter)\n",
    "gravity.particles.add_particles(bodies-moon)\n",
    "channel = {\"from stars\": bodies.new_channel_to(gravity.particles),\n",
    "            \"to_stars\": gravity.particles.new_channel_to(bodies)}\n",
    "\n",
    "hydro = Fi(converter, mode=\"openmp\")\n",
    "hydro.parameters.use_hydro_flag = True\n",
    "hydro.parameters.radiation_flag = False\n",
    "hydro.parameters.gamma = 1\n",
    "hydro.parameters.isothermal_flag = True\n",
    "hydro.parameters.integrate_entropy_flag = False\n",
    "hydro.parameters.timestep = 0.01*Pinner \n",
    "hydro.parameters.verbosity = 0\n",
    "hydro.parameters.eps_is_h_flag = False  # h_smooth is constant\n",
    "eps = 0.1 | units.au\n",
    "hydro.parameters.gas_epsilon = eps\n",
    "hydro.parameters.sph_h_const = eps\n",
    "\n",
    "hydro.particles.add_particles(disk)\n",
    "hydro.dm_particles.add_particles(moon.as_set())\n",
    "channel.update({\"from_disk\": disk.new_channel_to(hydro.particles)})\n",
    "channel.update({\"to_disk\": hydro.particles.new_channel_to(disk)})\n",
    "channel.update({\"from_moon\": moon.new_channel_to(hydro.dm_particles)})\n",
    "channel.update({\"to_moon\": hydro.dm_particles.new_channel_to(moon)})\n",
    "\n",
    "bodies.add_particles(disk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.couple import bridge\n",
    "from amuse.ext.composition_methods import *\n",
    "gravhydro = bridge.Bridge(use_threading=False)  # , method=SPLIT_4TH_S_M4)\n",
    "gravhydro.add_system(gravity, (hydro,))\n",
    "gravhydro.add_system(hydro, (gravity,))\n",
    "gravhydro.timestep = 0.2 * Pinner"
   ]
  },
  {
   "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_hydro_bridge(gravity, hydro, gravhydro, bodies,\n",
    "                         t_end):\n",
    "\n",
    "    gravity_initial_total_energy = gravity.get_total_energy() + hydro.get_total_energy()\n",
    "    model_time = 0 | units.Myr\n",
    "    dt = 0.012 | units.yr  # 1.0*Pinner\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:\", \"ae=\", orbit_planet[2].in_(units.AU), orbit_planet[3])\n",
    "        print(\"Moon:\", \"ae=\", orbit_moon[2].in_(units.AU), orbit_moon[3])\n",
    "        \n",
    "        dE_gravity = gravity_initial_total_energy / (\n",
    "            gravity.get_total_energy() + hydro.get_total_energy()\n",
    "        )\n",
    "        print(\"Time:\", model_time.in_(units.day), \\\n",
    "              \"dE=\", dE_gravity)  # , dE_hydro\n",
    "\n",
    "        gravhydro.evolve_model(model_time)\n",
    "        channel[\"to_stars\"].copy()\n",
    "        channel[\"to_disk\"].copy()\n",
    "        channel[\"to_moon\"].copy()\n",
    "        \n",
    "        print(\"S=\", bodies[:3])\n",
    "        print(\"g=\", gravity.particles)\n",
    "        print(gravity.particles.y.in_(units.au), moon.y.in_(units.au))\n",
    "\n",
    "    gravity.stop()\n",
    "    hydro.stop()\n",
    "\n",
    "t_end = 1.0 | units.yr\n",
    "gravity_hydro_bridge(gravity, hydro, gravhydro, \n",
    "                     bodies, t_end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(moon, disk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have adopted the star-planet-moon system with disk. Instead of hyrarchical_bridge.py where the disk was integrated using a direct N-body code, the disk is now a gaseous.\n",
    "\n",
    "\n",
    "Assignments and questions:\n",
    "---------------\n",
    "\n",
    "### Assignment 1:\n",
    "\n",
    "\n",
    "### Question 1:\n"
   ]
  },
  {
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}