1{
2 "cells": [
3  {
4   "cell_type": "markdown",
5   "id": "shaped-chicken",
6   "metadata": {},
7   "source": [
8    "Customising the adaptive integrator\n",
9    "===================================\n",
10    "\n",
11    "In the previous section we showed a few\n",
12    "usage examples of the ``taylor_adaptive`` class using the default\n",
13    "options. Here, we will show how the behaviour of the integrator\n",
14    "can be customised in a variety of ways.\n",
15    "\n",
16    "Error tolerance\n",
17    "---------------\n",
18    "\n",
19    "As we mentioned earlier, by default the ``taylor_adaptive`` class\n",
20    "uses an error tolerance equal to the machine\n",
21    "epsilon of the floating-point type in use. E.g., when using double precision,\n",
22    "the tolerance is set to $\\sim 2.2\\times 10^{-16}$.\n",
23    "\n",
24    "The tolerance value is used by the ``taylor_adaptive``\n",
25    "class to control the error arising from truncating the (infinite)\n",
26    "Taylor series representing the solution of the ODE system.\n",
27    "In other words, ``taylor_adaptive`` strives to ensure that the\n",
28    "magnitude of the remainders of the Taylor series is\n",
29    "not greater than the tolerance,\n",
30    "either in an absolute or relative sense. Absolute error control mode\n",
31    "is activated when all elements of the state vector have a magnitude\n",
32    "less than 1, while relative error control mode is activated when at least one\n",
33    "element of the state vector has a magnitude greater than 1.\n",
34    "\n",
35    "In order to specify a non-default tolerance, the keyword argument\n",
36    "``tol`` can be used when constructing an integrator object:"
37   ]
38  },
39  {
40   "cell_type": "code",
41   "execution_count": 1,
42   "id": "informational-liquid",
43   "metadata": {},
44   "outputs": [
45    {
46     "data": {
47      "text/plain": [
48       "Taylor order: 12\n",
49       "Dimension   : 2\n",
50       "Time        : 0.0000000000000000\n",
51       "State       : [0.050000000000000003, 0.025000000000000001]"
52      ]
53     },
54     "execution_count": 1,
55     "metadata": {},
56     "output_type": "execute_result"
57    }
58   ],
59   "source": [
60    "import heyoka as hy\n",
61    "\n",
62    "# Create the symbolic variables x and v.\n",
63    "x, v = hy.make_vars(\"x\", \"v\")\n",
64    "\n",
65    "# Create the integrator object.\n",
66    "ta = hy.taylor_adaptive(\n",
67    "                        # Definition of the ODE system:\n",
68    "                        # x' = v\n",
69    "                        # v' = -9.8 * sin(x)\n",
70    "                        sys = [(x, v),\n",
71    "                         (v, -9.8 * hy.sin(x))],\n",
72    "                        # Initial conditions for x and v.\n",
73    "                        state = [0.05, 0.025],\n",
74    "                        # Set the tolerance to 1e-9\n",
75    "                        tol = 1e-9)\n",
76    "\n",
77    "ta"
78   ]
79  },
80  {
81   "cell_type": "markdown",
82   "id": "confirmed-majority",
83   "metadata": {},
84   "source": [
85    "The optimal Taylor order for a tolerance of $10^{-9}$\n",
86    "is now 12 (instead of 20 for a tolerance of $\\sim 2.2\\times 10^{-16}$).\n",
87    "\n",
88    "Integrating the system back and forth shows how the accuracy of the\n",
89    "integration is reduced with respect to the default tolerance value:"
90   ]
91  },
92  {
93   "cell_type": "code",
94   "execution_count": 2,
95   "id": "technological-holiday",
96   "metadata": {},
97   "outputs": [
98    {
99     "data": {
100      "text/plain": [
101       "Taylor order: 12\n",
102       "Dimension   : 2\n",
103       "Time        : 0.0000000000000000\n",
104       "State       : [0.050000000001312876, 0.024999999997558964]"
105      ]
106     },
107     "execution_count": 2,
108     "metadata": {},
109     "output_type": "execute_result"
110    }
111   ],
112   "source": [
113    "# Integrate forth to t = 10 and then back to t = 0.\n",
114    "ta.propagate_until(t = 10.)\n",
115    "ta.propagate_until(t = 0.)\n",
116    "\n",
117    "ta"
118   ]
119  },
120  {
121   "cell_type": "markdown",
122   "id": "national-investing",
123   "metadata": {},
124   "source": [
125    "Compact mode\n",
126    "------------\n",
127    "\n",
128    "By default, the just-in-time compilation process of heyoka.py\n",
129    "aims at maximising runtime performance over everything else.\n",
130    "In practice, this means that heyoka.py generates a timestepper\n",
131    "function in which there are no branches and where all loops\n",
132    "have been fully unrolled.\n",
133    "\n",
134    "This approach leads to highly optimised timestepper functions,\n",
135    "but, on the other hand, it can result in long compilation times\n",
136    "and high memory usage for large ODE systems. Thus, heyoka.py provides\n",
137    "also a *compact mode* option in which code generation employs\n",
138    "more traditional programming idioms that greatly reduce compilation time\n",
139    "and memory usage. Compact mode results in a performance degradation\n",
140    "of $\\lesssim 2\\times$ with respect to the default code generation\n",
141    "mode, but it renders heyoka.py usable with ODE systems consisting\n",
142    "of thousands of terms.\n",
143    "\n",
144    "Let's try to quantify the performance difference in a concrete case.\n",
145    "In this example, we first construct the ODE system corresponding\n",
146    "to an N-body problem with 6 particles via the ``make_nbody_sys()``\n",
147    "utility function:"
148   ]
149  },
150  {
151   "cell_type": "code",
152   "execution_count": 3,
153   "id": "employed-bruce",
154   "metadata": {},
155   "outputs": [],
156   "source": [
157    "# Create an nbody system with 6 particles.\n",
158    "sys = hy.make_nbody_sys(n = 6)"
159   ]
160  },
161  {
162   "cell_type": "markdown",
163   "id": "seasonal-three",
164   "metadata": {},
165   "source": [
166    "Next, we create an initial state vector for our system.\n",
167    "The contents of the vector do not matter at this stage:"
168   ]
169  },
170  {
171   "cell_type": "code",
172   "execution_count": 4,
173   "id": "suburban-cannon",
174   "metadata": {},
175   "outputs": [],
176   "source": [
177    "# Create an initial state vector (6 values per body).\n",
178    "import numpy as np\n",
179    "sv = np.zeros(36)"
180   ]
181  },
182  {
183   "cell_type": "markdown",
184   "id": "western-metallic",
185   "metadata": {},
186   "source": [
187    "Next, we time the creation of an integrator object in default\n",
188    "code generation mode:"
189   ]
190  },
191  {
192   "cell_type": "code",
193   "execution_count": 5,
194   "id": "super-jimmy",
195   "metadata": {},
196   "outputs": [
197    {
198     "name": "stdout",
199     "output_type": "stream",
200     "text": [
201      "CPU times: user 2.93 s, sys: 17.7 ms, total: 2.95 s\n",
202      "Wall time: 2.96 s\n"
203     ]
204    }
205   ],
206   "source": [
207    "%time ta_default = hy.taylor_adaptive(sys, sv)"
208   ]
209  },
210  {
211   "cell_type": "markdown",
212   "id": "shared-exhaust",
213   "metadata": {},
214   "source": [
215    "Finally, we time the creation of the same integrator object\n",
216    "in compact mode (which can be activated via the ``compact_mode``\n",
217    "keyword argument):"
218   ]
219  },
220  {
221   "cell_type": "code",
222   "execution_count": 6,
223   "id": "elder-metadata",
224   "metadata": {},
225   "outputs": [
226    {
227     "name": "stdout",
228     "output_type": "stream",
229     "text": [
230      "CPU times: user 206 ms, sys: 580 µs, total: 207 ms\n",
231      "Wall time: 210 ms\n"
232     ]
233    }
234   ],
235   "source": [
236    "%time ta_default = hy.taylor_adaptive(sys, sv, compact_mode = True)"
237   ]
238  },
239  {
240   "cell_type": "markdown",
241   "id": "funny-legislature",
242   "metadata": {},
243   "source": [
244    "That is, in this specific example compact mode is more than 10 times\n",
245    "faster than the default\n",
246    "code generation mode when it comes to the construction of the integrator\n",
247    "object. For larger ODE systems, the gap will be even wider.\n",
248    "\n",
249    "High-accuracy mode\n",
250    "------------------\n",
251    "\n",
252    "For long-term integrations at very low error tolerances, heyoka.py offers\n",
253    "an opt-in *high-accuracy* mode. In high-accuracy mode, heyoka.py\n",
254    "employs techniques that minimise the numerical errors arising from\n",
255    "the use of finite-precision floating-point numbers, at the cost\n",
256    "of a slight runtime performance degradation.\n",
257    "\n",
258    "Currently, high-accuracy mode changes the way heyoka,py evaluates\n",
259    "the Taylor polynomials used to update the state of the system\n",
260    "at the end of an integration timestep. Specifically, polynomial evaluation\n",
261    "via Horner's rule is replaced by\n",
262    "[compensated summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm),\n",
263    "which prevents catastrophic cancellation issues and ultimately helps maintaining\n",
264    "machine precision over very long integrations.\n",
265    "\n",
266    "High-accuracy mode can be enabled via the ``high_accuracy`` keyword\n",
267    "argument."
268   ]
269  }
270 ],
271 "metadata": {
272  "kernelspec": {
273   "display_name": "Python 3",
274   "language": "python",
275   "name": "python3"
276  },
277  "language_info": {
278   "codemirror_mode": {
279    "name": "ipython",
280    "version": 3
281   },
282   "file_extension": ".py",
283   "mimetype": "text/x-python",
284   "name": "python",
285   "nbconvert_exporter": "python",
286   "pygments_lexer": "ipython3",
287   "version": "3.8.6"
288  }
289 },
290 "nbformat": 4,
291 "nbformat_minor": 5
292}
293