diff --git a/ipython_examples/AdvWHFast.ipynb b/ipython_examples/AdvWHFast.ipynb index de3907812..3d423c865 100644 --- a/ipython_examples/AdvWHFast.ipynb +++ b/ipython_examples/AdvWHFast.ipynb @@ -2,7 +2,10 @@ "cells": [ { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "# Advanced settings for WHFast: Extra speed, accuracy, and additional forces\n", "\n", @@ -22,7 +25,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "**The Wisdom-Holman algorithm**\n", "\n", @@ -39,7 +45,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "**Combining Kepler steps and synchronizing**\n", "\n", @@ -55,7 +64,7 @@ "\n", "**Conversions between Jacobi and Inertial Coordinates**\n", "\n", - "It turns out that the most convenient coordinate system to work in for performing the Kepler steps is Jacobi coordinates (see, e.g., 9.5.4 of Murray & Dermott). WHFast therefore works in Jacobi coordinates, converting to inertial coordinates when it needs to (e.g. for output, and for doing the direct gravity calculation in the interaction step, which is most easily done in inertial coordinates).\n", + "It turns out that the most convenient coordinate system to work in for performing the Kepler steps is often Jacobi coordinates (see, e.g., 9.5.4 of Murray & Dermott). WHFast therefore works in Jacobi coordinates by default, converting to inertial coordinates when it needs to (e.g. for output, and for doing the direct gravity calculation in the interaction step, which is most easily done in inertial coordinates).\n", "\n", "One feature of WHFast is that it works in whatever inertial coordinate system you choose for your initial conditions. This means that whatever happens behind the scenes, the user always gets the particles' inertial coordinates at the front end. At the beginning of every timestep, WHFast therefore has to somehow obtain the Jacobi coordinates. The straightforward thing would be to convert from the inertial coordinates to Jacobi coordinates every timestep, but these conversions slow things down, and they represent extra operations that grow the round-off error.\n", "\n", @@ -68,9 +77,11 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "metadata": { - "collapsed": false + "collapsed": false, + "deletable": true, + "editable": true }, "outputs": [], "source": [ @@ -88,7 +99,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "By default WHFast synchronizes and recalculates the Jacobi coordinates from the inertial ones every timestep. This guarantees that the user always gets physical particle states for output, and ensures reliable output if the user decides to, e.g., grow the particles' masses between timesteps. \n", "\n", @@ -97,9 +111,11 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": { - "collapsed": false + "collapsed": false, + "deletable": true, + "editable": true }, "outputs": [], "source": [ @@ -109,16 +125,21 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "Now it becomes the user's responsibility to appropriately synchronize and recalculate jacobi coordinates when needed. You can tell WHFast to recalculate Jacobi coordinates for a given timestep (say after you change a particle's mass) with the `sim.ri_whfast.recalculate_jacobi_this_timestep` flag. After it recalculates Jacobi coordinates, WHFast will reset this flag to zero, so you just set it each time you mess with the particles." ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "metadata": { - "collapsed": false + "collapsed": false, + "deletable": true, + "editable": true }, "outputs": [ { @@ -127,29 +148,29 @@ "text": [ "safe_mode = 1\n", "---------------------------------\n", - "REBOUND version: \t2.18.7\n", - "REBOUND built on: \tJun 20 2016 19:21:51\n", + "REBOUND version: \t3.4.0\n", + "REBOUND built on: \tMay 31 2017 11:53:50\n", "Number of particles: \t2\n", "Selected integrator: \twhfast\n", - "Simulation time: \t628318.530718\n", + "Simulation time: \t6.2831853071795858e+05\n", "Current timestep: \t0.200000\n", "---------------------------------\n", - "\n", - "\n", + "\n", + "\n", "---------------------------------\n", - "Safe integration took 1.92633914948 seconds\n", + "Safe integration took 1.4043679237365723 seconds\n", "---------------------------------\n", - "REBOUND version: \t2.18.7\n", - "REBOUND built on: \tJun 20 2016 19:21:51\n", + "REBOUND version: \t3.4.0\n", + "REBOUND built on: \tMay 31 2017 11:53:50\n", "Number of particles: \t2\n", "Selected integrator: \twhfast\n", - "Simulation time: \t628318.530718\n", + "Simulation time: \t6.2831853071795858e+05\n", "Current timestep: \t0.200000\n", "---------------------------------\n", - "\n", - "\n", + "\n", + "\n", "---------------------------------\n", - "Manual integration took 1.18098711967 seconds\n" + "Manual integration took 0.8836901187896729 seconds\n" ] } ], @@ -175,7 +196,9 @@ { "cell_type": "markdown", "metadata": { - "collapsed": true + "collapsed": true, + "deletable": true, + "editable": true }, "source": [ "In our test case with a single planet, there is effectively no interaction step, and by combining Kepler steps we get almost the full factor of 2 speedup we expect. Because Kepler steps are expensive (by virtue of having to solve the transcendental Kepler equation), this will always be an important performance boost for few-planet cases.\n", @@ -185,9 +208,11 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": { - "collapsed": true + "collapsed": true, + "deletable": true, + "editable": true }, "outputs": [], "source": [ @@ -196,7 +221,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "REBOUND will synchronize every timestep even if you set `sim.ri_whfast.safe_mode = 0` and never explicitly call `sim.integrator_synchronize()`." ] @@ -204,7 +232,9 @@ { "cell_type": "markdown", "metadata": { - "collapsed": true + "collapsed": true, + "deletable": true, + "editable": true }, "source": [ "**Modifying particles/forces**\n", @@ -214,9 +244,11 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": { - "collapsed": true + "collapsed": true, + "deletable": true, + "editable": true }, "outputs": [], "source": [ @@ -232,7 +264,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "Here, because we grow the mass of the planet every timestep, we have to recalculate Jacobi coordinates every timestep (since they depend on the masses of the particles). We therefore manually set the flag to recalculate them the next timestep every time we make a change. Here we would actually get the same result if we just left `sim.ri_whfast.safe_mode = 1`, since when recalculating Jacobi coordinates, WHFast automatically has to synchronize in order to get real positions and velocities for the planets. In this case WHFast is therefore synchronizing and recalculating Jacobi coordinates every timestep.\n", "\n", @@ -241,9 +276,11 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": { - "collapsed": false + "collapsed": false, + "deletable": true, + "editable": true }, "outputs": [], "source": [ @@ -259,14 +296,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "This would not give accurate results, because the `sim.particles[1].vx` we access after `sim.step()` isn't a physical velocity (it's missing a half-Kepler step). It's basically at an intermediate point in the calculation. In order to make this work, one would call `sim.integrator_synchronize()` between `sim.step()` and accessing `sim.particles[1].vx`, to ensure the velocity is physical." ] }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "**Symplectic correctors**\n", "\n", @@ -275,9 +318,11 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, "metadata": { - "collapsed": true + "collapsed": true, + "deletable": true, + "editable": true }, "outputs": [], "source": [ @@ -286,7 +331,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "By default, WHFast does not use correctors, i.e., sim.integrator_whfast_corrector = 0. This is because the default is also to synchronize every timestep. An Nth order corrector does N-1 Kepler steps of various sizes, so an 11th order corrector done every timestep would increase the number of Kepler steps by an order of magnitude, making WHFast unacceptably slow. So keep in mind that if you're doing modifications that require recalculating jacobi coordinates or synchronizing every timestep, you should turn off symplectic correctors (the default) unless you really need the accuracy." ] @@ -294,6 +342,38 @@ { "cell_type": "markdown", "metadata": {}, + "source": [ + "**Changing the internal coordinate system**\n", + "\n", + "WHFast by default uses Jacobi coordinates internally. This works well for planetary systems which are stable and orbits are not crossing. However, in some cases a different coordinate system might perform better. WHFast also support so called democratic heliocentric coordinates and the so called WHDS coordinates. For mor information on these coordinates systems [see Hernandez and Dehnen (2016)](https://arxiv.org/abs/1612.05329). To select a different coordinate system, use the following syntax:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "sim.ri_whfast.coordinates = 'jacobi' #default\n", + "sim.ri_whfast.coordinates = 'democraticheliocentric' \n", + "sim.ri_whfast.coordinates = 'whds' " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that symplectic corrector are only compatible with Jacobi coordinates because both democratic heliocentric and WHDS include a so called jump step." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "**Warning messages**\n", "\n", @@ -304,7 +384,9 @@ "cell_type": "code", "execution_count": 8, "metadata": { - "collapsed": false + "collapsed": false, + "deletable": true, + "editable": true }, "outputs": [ { @@ -326,7 +408,9 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true + "collapsed": true, + "deletable": true, + "editable": true }, "outputs": [], "source": [] @@ -348,7 +432,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.5.1" + "version": "3.5.2" } }, "nbformat": 4,