Add IAS15 high-precision integrator (#643)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Jo Bovy <bovy@astro.utoronto.ca>

Author: 3 people

HISTORY.txt

@@ -1,6 +1,8 @@
 v1.10.2 (expected around 2025-03-01)
 ====================
 
+ - Implemented the IAS15 integrator of Rein & Spiegel 2015.
+
  - Combine the drift calculations in the Python leapfrog integrator (#365).
 
  - Move the checks for non-axisymmetric and dissipative potentials from internal to

doc/source/orbit.rst

@@ -1014,9 +1014,11 @@ the ``orbit.integrate`` method. Currently available integrators are
 * rk6_c
 * dopr54_c
 * dop853_c
+* ias15_c
 
-which are Runge-Kutta and Dormand-Prince methods. There are also a
-number of symplectic integrators available
+which are Runge-Kutta, Dormand-Prince methods and the IAS15
+integrator in `Rein & Spiegel (2014) <https://doi.org/10.1093/mnras/stu2164>`_.
+There are also a number of symplectic integrators available
 
 * leapfrog_c
 * symplec4_c
@@ -1045,6 +1047,9 @@ which are speedy and reliable. For example, compare
 >>> timeit(o.integrate(ts,mp,method='dop853_c'))
 # 4.65 ms ± 86.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
 
+The ``ias15_c`` method uses adaptive timestepping under the hood, and
+can be used in cases where very high precision is required during orbit integration.
+
 If the C extensions are unavailable for some reason, I recommend using
 the ``odeint`` pure-Python integrator, as it is the fastest. Using the
 same example as above

galpy/orbit/Orbits.py

@@ -1287,13 +1287,15 @@ def check_integrator(method, no_symplec=False):
             "rk6_c",
             "dopr54_c",
             "dop853_c",
+            "ias15_c",
         ]
         if no_symplec:
             symplec_methods = [
                 "leapfrog",
                 "leapfrog_c",
                 "symplec4_c",
                 "symplec6_c",
+                "ias15_c",  # practically speaking, ias15 has the same limitations as symplectic integrators in galpy
             ]
             [valid_methods.remove(symplec_method) for symplec_method in symplec_methods]
         if method.lower() not in valid_methods:
@@ -1385,6 +1387,7 @@ def integrate(
           -  'dopr54_c' for a 5-4 Dormand-Prince integrator in C
           -  'dop853' for a 8-5-3 Dormand-Prince integrator in Python
           -  'dop853_c' for a 8-5-3 Dormand-Prince integrator in C
+          -  'ias15_c' for an adaptive 15th order integrator using Gauß-Radau quadrature (see IAS15 paper) in C
 
         - 2018-10-13 - Written as parallel_map applied to regular Orbit integration - Mathew Bub (UofT)
         - 2018-12-26 - Written to use OpenMP C implementation - Bovy (UofT)
@@ -5153,6 +5156,7 @@ def bruteSOS(
           -  'dopr54_c' for a 5-4 Dormand-Prince integrator in C
           -  'dop853' for a 8-5-3 Dormand-Prince integrator in Python
           -  'dop853_c' for a 8-5-3 Dormand-Prince integrator in C
+          -  'ias15_c' for an adaptive 15th order integrator using Gauß-Radau quadrature (see IAS15 paper) in C
 
         - 2023-05-31 - Written - Bovy (UofT)
 
@@ -5895,6 +5899,7 @@ def plotBruteSOS(
           -  'dopr54_c' for a 5-4 Dormand-Prince integrator in C
           -  'dop853' for a 8-5-3 Dormand-Prince integrator in Python
           -  'dop853_c' for a 8-5-3 Dormand-Prince integrator in C
+          -  'ias15_c' for an adaptive 15th order integrator using Gauß-Radau quadrature (see IAS15 paper) in C
 
         - 2023-05-31 - Written - Bovy (UofT)
 

galpy/orbit/integratePlanarOrbit.py

@@ -618,6 +618,8 @@ def _parse_integrator(int_method):
         int_method_c = 5
     elif int_method.lower() == "dop853_c":
         int_method_c = 6
+    elif int_method.lower() == "ias15_c":
+        int_method_c = 7
     else:
         int_method_c = 0
     return int_method_c

galpy/orbit/orbit_c_ext/integrateFullOrbit.c

@@ -14,6 +14,7 @@
 #include <bovy_symplecticode.h>
 #include <leung_dop853.h>
 #include <bovy_rk.h>
+#include <wez_ias15.h>
 #include <integrateFullOrbit.h>
 //Potentials
 #include <galpy_potentials.h>
@@ -717,6 +718,11 @@ EXPORT void integrateFullOrbit(int nobj,
     odeint_deriv_func= &evalRectDeriv;
     dim= 6;
     break;
+  case 7: //ias15
+    odeint_func= &wez_ias15;
+    odeint_deriv_func= &evalRectForce;
+    dim= 3;
+    break;
   }
 #pragma omp parallel for schedule(dynamic,ORBITS_CHUNKSIZE) private(ii,jj) num_threads(max_threads)
   for (ii=0; ii < nobj; ii++) {
@@ -883,6 +889,11 @@ void integrateOrbit_dxdv(double *yo,
     odeint_deriv_func= &evalRectDeriv_dxdv;
     dim= 12;
     break;
+  case 7: //ias15
+    odeint_func= &wez_ias15;
+    odeint_deriv_func= &evalRectForce;
+    dim= 6;
+    break;
   }
   odeint_func(odeint_deriv_func,dim,yo,nt,-9999.99,t,npot,potentialArgs,
 	      rtol,atol,result,err);

galpy/orbit/orbit_c_ext/integrateLinearOrbit.c

@@ -7,6 +7,7 @@
 #include <math.h>
 #include <bovy_symplecticode.h>
 #include <bovy_rk.h>
+#include <wez_ias15.h>
 #include <leung_dop853.h>
 #include <integrateFullOrbit.h>
 //Potentials
@@ -192,6 +193,11 @@ EXPORT void integrateLinearOrbit(int nobj,
     odeint_deriv_func= &evalLinearDeriv;
     dim= 2;
     break;
+  case 7: //ias15
+    odeint_func= &wez_ias15;
+    odeint_deriv_func= &evalLinearForce;
+    dim= 1;
+    break;
   }
 #pragma omp parallel for schedule(dynamic,ORBITS_CHUNKSIZE) private(ii) num_threads(max_threads)
   for (ii=0; ii < nobj; ii++) {

galpy/orbit/orbit_c_ext/integratePlanarOrbit.c

@@ -8,6 +8,7 @@
 #include <bovy_coords.h>
 #include <bovy_symplecticode.h>
 #include <bovy_rk.h>
+#include <wez_ias15.h>
 #include <leung_dop853.h>
 #include <integrateFullOrbit.h>
 //Potentials
@@ -674,6 +675,11 @@ EXPORT void integratePlanarOrbit(int nobj,
     odeint_deriv_func= &evalPlanarRectDeriv;
     dim= 4;
     break;
+  case 7: //ias15
+    odeint_func= &wez_ias15;
+    odeint_deriv_func= &evalPlanarRectForce;
+    dim= 2;
+    break;
   }
 #pragma omp parallel for schedule(dynamic,ORBITS_CHUNKSIZE) private(ii,jj) num_threads(max_threads)
   for (ii=0; ii < nobj; ii++) {