added new adaptive timestep scheme from REBOUND (PRS24)

rahil-makadia · Mar 13, 2024 · b5be3df · b5be3df
1 parent 0cc9511
commit b5be3df
Show file tree

Hide file tree

Showing 9 changed files with 73 additions and 34 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -16,7 +16,7 @@ project(grss VERSION ${ver})
 set(CMAKE_CXX_STANDARD 11)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 add_compile_options(-std=c++11 -O3 -fPIC -fopenmp) # operational flags
-# add_compile_options(-std=c++11 -DLONGDOUBLE -g3 -fPIC -fopenmp -Werror -Wall -Wextra -pedantic) # debugging flags
+# add_compile_options(-std=c++11 -g3 -fPIC -fopenmp -Werror -Wall -Wextra -pedantic) # debugging flags
 
 # Set header file directories
 include_directories(${CMAKE_SOURCE_DIR}/include)

diff --git a/grss/version.txt b/grss/version.txt
@@ -1 +1 @@
-3.4.2
+3.5.0
diff --git a/include/interpolate.h b/include/interpolate.h
@@ -3,7 +3,7 @@
 
 #include "force.h"
 
-inline void comp_sum(real num, real *sum, real *compCoeff) {
+static inline void comp_sum(real num, real *sum, real *compCoeff) {
     const real y = num - *compCoeff;
     const real t = *sum + y;
     *compCoeff = (t - *sum) - y;

diff --git a/joss/joss_paper.bib b/joss/joss_paper.bib
@@ -89,7 +89,7 @@ @article{Holman2023
 }
 
 @article{Schwamb2023,
-author    = {Schwamb, Megan E. and Jones, R. Lynne and Yoachim, Peter and Volk, Kathryn and Dorsey, Rosemary C. and Opitom, Cyrielle and Greenstreet, Sarah and Lister, Tim and Snodgrass, Colin and Bolin, Bryce T. and Inno, Laura and Bannister, Michele T. and Eggl, Siegfried and Solontoi, Michael and Kelley, Michael S. P. and Jurić, Mario and Lin 林, Hsing Wen 省文 and Ragozzine, Darin and Bernardinelli, Pedro H. and Chesley, Steven R. and Daylan, Tansu and Ďurech, Josef and Fraser, Wesley C. and Granvik, Mikael and Knight, Matthew M. and Lisse, Carey M. and Malhotra, Renu and Oldroyd, William J. and Thirouin, Audrey and Ye 叶, Quanzhi 泉志},
+author    = {Schwamb, Megan E. and Jones, R. Lynne and Yoachim, Peter and Volk, Kathryn and Dorsey, Rosemary C. and Opitom, Cyrielle and Greenstreet, Sarah and Lister, Tim and Snodgrass, Colin and Bolin, Bryce T. and Inno, Laura and Bannister, Michele T. and Eggl, Siegfried and Solontoi, Michael and Kelley, Michael S. P. and Jurić, Mario and Lin, Hsing Wen and Ragozzine, Darin and Bernardinelli, Pedro H. and Chesley, Steven R. and Daylan, Tansu and Ďurech, Josef and Fraser, Wesley C. and Granvik, Mikael and Knight, Matthew M. and Lisse, Carey M. and Malhotra, Renu and Oldroyd, William J. and Thirouin, Audrey and Ye, Quanzhi},
 journal   = {The Astrophysical Journal Supplement Series},
 title     = {Tuning the Legacy Survey of Space and Time (LSST) Observing Strategy for Solar System Science},
 year      = {2023},

diff --git a/src/gr15.cpp b/src/gr15.cpp
@@ -217,6 +217,54 @@ static real root7(real num){
     return fac*root;
 }
 
+// static real get_adaptive_timestep_old(propSimulation *propSim, const real &dt,
+//                                       const real &accIntegArr7Max,
+//                                       const std::vector<std::vector<real>> &b) {
+//     real b6Max;
+//     vabs_max(b[6], b6Max);
+//     real relError = b6Max / accIntegArr7Max;
+//     real dtReq;
+//     if (std::isnormal(relError)) {
+//         dtReq = root7(propSim->integParams.tolInteg/relError)*dt;
+//     } else {
+//         dtReq = dt/propSim->integParams.dtChangeFactor;
+//     }
+//     return dtReq;
+// }
+
+static real get_adaptive_timestep(propSimulation *propSim, const real &dt,
+                                  const std::vector<real> &accInteg0,
+                                  const std::vector<std::vector<real>> &b) {
+    real minTimescale2 = 1e300L;
+    size_t startb = 0;
+    for (size_t i = 0; i < propSim->integParams.nInteg; i++){
+        real y2 = 0.0L;
+        real y3 = 0.0L;
+        real y4 = 0.0L;
+        for (size_t j = 0; j < 3; j++){
+            const size_t k = startb+j;
+            real temp = accInteg0[k] + b[0][k] + b[1][k] + b[2][k] + b[3][k] + b[4][k] + b[5][k] + b[6][k];
+            y2 += temp*temp;
+            temp = b[0][k] + 2.0L*b[1][k] + 3.0L*b[2][k] + 4.0L*b[3][k] + 5.0L*b[4][k] + 6.0L*b[5][k] + 7.0L*b[6][k];
+            y3 += temp*temp;
+            temp = 2.0L*b[1][k] + 6.0L*b[2][k] + 12.0L*b[3][k] + 20.0L*b[4][k] + 30.0L*b[5][k] + 42.0L*b[6][k];
+            y4 += temp*temp;
+            real timescale2 = 2.0L*y2/(y3+sqrt(y4*y2));
+            if (std::isnormal(timescale2) && timescale2 < minTimescale2){
+                minTimescale2 = timescale2;
+            }
+        }
+    }
+    real dtReq;
+    if (std::isnormal(minTimescale2)){
+        // Numerical factor below is from REBOUND
+        dtReq = sqrt(minTimescale2) * root7(propSim->integParams.tolInteg*5040.0) * dt;
+    }else{
+        dtReq = dt/propSim->integParams.dtChangeFactor;
+    }
+    return dtReq;
+}
+
 void gr15(propSimulation *propSim) {
     if (!std::isfinite(propSim->t)) {
         throw std::runtime_error("t is not finite");
@@ -249,8 +297,7 @@ void gr15(propSimulation *propSim) {
     real *b6Tilde = new real[dim];
     memset(b6Tilde, 0.0, dim * sizeof(real));
     real b6TildeMax, accIntegArr7Max;
-    real b6Max, accIntegNextMax;
-    real relError, dtReq;
+    real dtReq;
     real tNextEvent = propSim->integParams.tf;
     static size_t nextEventIdx = 0;
     if (t == propSim->integParams.t0) {
@@ -309,27 +356,20 @@ void gr15(propSimulation *propSim) {
                 PCerr = b6TildeMax / accIntegArr7Max;
                 bOld = b;
             }
-            vabs_max(b[6], b6Max);
-            vabs_max(accIntegArr[7], accIntegNextMax);
-            relError = b6Max / accIntegNextMax;
             if (propSim->integParams.adaptiveTimestep) {
-                if (std::isnormal(relError)) {
-                    dtReq = root7(propSim->integParams.tolInteg/relError)*dt;
-                } else {
-                    dtReq = dt/propSim->integParams.dtChangeFactor;
-                }
+                // dtReq = get_adaptive_timestep_old(propSim, dt, accIntegArr7Max, b);
+                dtReq = get_adaptive_timestep(propSim, dt, accInteg0, b);
             } else {
                 dtReq = dt;
             }
             if (fabs(dtReq) < propSim->integParams.dtMin) {
                 dtReq = copysign(propSim->integParams.dtMin, dtReq);
             }
             if (fabs(dtReq/dt) < propSim->integParams.dtChangeFactor) {
-                if (propSim->integParams.timestepCounter > 1) {
-                    refine_b(b, e, dtReq / dt, dim);
-                }
+                // std::cout << "Restarting next while loop at time t = " << t
+                //           << " with dt = " << dtReq << " instead of " << dt
+                //           << std::endl;
                 dt = dtReq;
-                std::cout << "Restarting next while loop at time t = " << t << std::endl;
                 continue;
             }
             // accept step

diff --git a/tests/python/prop/apophis.ipynb b/tests/python/prop/apophis.ipynb
diff --git a/tests/python/prop/didymos.ipynb b/tests/python/prop/didymos.ipynb
@@ -168,7 +168,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.9"
+   "version": "3.10.13"
   },
   "orig_nbformat": 4
  },

diff --git a/tests/python/prop/icarus.ipynb b/tests/python/prop/icarus.ipynb
@@ -167,7 +167,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.9"
+   "version": "3.10.13"
   },
   "orig_nbformat": 4
  },

diff --git a/tests/python/prop/moshup.ipynb b/tests/python/prop/moshup.ipynb