From dec6d786906d3064b9feca9c696114f88e124762 Mon Sep 17 00:00:00 2001 From: Hanno Rein Date: Mon, 6 May 2024 08:12:40 -0400 Subject: [PATCH] switching function --- docs/index.md | 1 + docs/integrators.md | 14 ++++++++++---- rebound/integrators/trace.py | 2 +- src/integrator_trace.c | 8 ++++---- src/rebound.h | 2 +- 5 files changed, 17 insertions(+), 10 deletions(-) diff --git a/docs/index.md b/docs/index.md index 845ff3406..5f7dec029 100644 --- a/docs/index.md +++ b/docs/index.md @@ -110,3 +110,4 @@ REBOUND is free software: you can redistribute it and/or modify it under the ter REBOUND is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with REBOUND. If not, see . + diff --git a/docs/integrators.md b/docs/integrators.md index 06c3f1ab0..bde2bd59a 100644 --- a/docs/integrators.md +++ b/docs/integrators.md @@ -404,7 +404,7 @@ The `reb_integrator_mercurius` structure contains the configuration and data str TRACE is a hybrid time-reversible integrator, based on the algorithm described in [Hernandez & Dehnen 2023](https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.4639H/abstract). It uses WHFast for long term integrations but switches time-reversibly to BS or IAS15 for all close encounters. TRACE is appropriate for systems with a dominant central mass that will occasionally have close encounters. -The TRACE implementation is described in UPDATE LINK [Lu et al in prep](https://ui.adsabs.harvard.edu/abs/2023MNRAS.522.4639H/abstract) +A paper describing the TRACE implementation is in preparation. The following code enables TRACE and sets the critical radius to 4 Hill radii @@ -429,7 +429,8 @@ The `reb_integrator_trace` structure contains the configuration and data structu If NULL (the default), the default switching function will be used. The arguments `i` and `j` are the indices of the two particles considered. The return value is a scalar. - If this function returns a negative value, it means a close encounter has been flagged. + A negative return value means a close encounter has been flagged. + If the return values of both this function and the central switching function below are always positive, then the integrator effectively becomes the standard Wisdom-Holman integrator. - Default switching function @@ -458,14 +459,14 @@ The `reb_integrator_trace` structure contains the configuration and data structu If NULL (the default), the default switching function will be used. The argument `j` is the index of the non-central particle considered. The return value is a scalar. - If this function returns a negative value, it means a close encounter has been flagged. If the return values of both this function and the non-central switching function are always positive, then the integrator effectively becomes the standard Wisdom-Holman integrator. + A negative value means a close encounter has been flagged. - Default switching function This switching function checks if a body is close to its pericenter by considering a timescale derived from high-order derivatives of the particle's herliocentric position, inspired by [Pham, Rein, and Spiegel 2024](https://ui.adsabs.harvard.edu/abs/2024OJAp....7E...1P/abstract). ```c - double reb_integrator_trace_switch_peri_pham(const struct reb_simulation* const r, const unsigned int j); + double reb_integrator_trace_switch_peri_default(const struct reb_simulation* const r, const unsigned int j); ``` - Fdot switching function @@ -488,14 +489,19 @@ The `reb_integrator_trace` structure contains the configuration and data structu === "C" ```c struct reb_simulation* r = reb_create_simulation(); + r->ri_trace.S_peri = reb_integrator_trace_switch_peri_default; // default r->ri_trace.S_peri = reb_integrator_trace_switch_peri_fdot; r->ri_trace.S_peri = reb_integrator_trace_switch_peri_distance; + r->ri_trace.S_peri = reb_integrator_trace_switch_peri_none; // Turn off pericenter switching ``` === "Python" ```python sim = rebound.Simulation() + sim.ri_trace.S_peri = "default" # Following Pham et al 2024 + sim.ri_trace.S_peri = "fdot" sim.ri_trace.S_peri = "distance" + sim.ri_trace.S_peri = "none" # Turn off pericenter switching ``` `double r_crit_hill` diff --git a/rebound/integrators/trace.py b/rebound/integrators/trace.py index 72ba658a2..6f979be7d 100644 --- a/rebound/integrators/trace.py +++ b/rebound/integrators/trace.py @@ -85,7 +85,7 @@ def S_peri(self): @S_peri.setter def S_peri(self, func): if func == "default": - self._S_peri = cast(clibrebound.reb_integrator_trace_switch_peri_pham,TRACECF) + self._S_peri = cast(clibrebound.reb_integrator_trace_switch_peri_default,TRACECF) elif func == "fdot": self._S_peri = cast(clibrebound.reb_integrator_trace_switch_peri_fdot,TRACECF) elif func == "distance": diff --git a/src/integrator_trace.c b/src/integrator_trace.c index b18650f20..cdf27dd60 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -115,8 +115,8 @@ int reb_integrator_trace_switch_peri_fdot(struct reb_simulation* const r, const return peff2 < pfdot*pfdot * r->dt*r->dt; } -int reb_integrator_trace_switch_peri_pham(struct reb_simulation* const r, const unsigned int j){ - // Many square roots, can this be fixed? +int reb_integrator_trace_switch_peri_default(struct reb_simulation* const r, const unsigned int j){ + // Following Pham et al (2024) const struct reb_integrator_trace* const ri_trace = &(r->ri_trace); double GM = r->G*r->particles[0].m; // Not sure if this is the right mass to use. @@ -544,7 +544,7 @@ void reb_integrator_trace_pre_ts_check(struct reb_simulation* const r){ const int N = r->N; const int Nactive = r->N_active==-1?r->N:r->N_active; int (*_switch) (struct reb_simulation* const r, const unsigned int i, const unsigned int j) = ri_trace->S ? ri_trace->S : reb_integrator_trace_switch_default; - int (*_switch_peri) (struct reb_simulation* const r, const unsigned int j) = ri_trace->S_peri ? ri_trace->S_peri : reb_integrator_trace_switch_peri_pham; + int (*_switch_peri) (struct reb_simulation* const r, const unsigned int j) = ri_trace->S_peri ? ri_trace->S_peri : reb_integrator_trace_switch_peri_default; // Clear encounter map for (unsigned int i=1; iN; i++){ @@ -621,7 +621,7 @@ double reb_integrator_trace_post_ts_check(struct reb_simulation* const r){ const int N = r->N; const int Nactive = r->N_active==-1?r->N:r->N_active; int (*_switch) (struct reb_simulation* const r, const unsigned int i, const unsigned int j) = ri_trace->S ? ri_trace->S : reb_integrator_trace_switch_default; - int (*_switch_peri) (struct reb_simulation* const r, const unsigned int j) = ri_trace->S_peri ? ri_trace->S_peri : reb_integrator_trace_switch_peri_pham; + int (*_switch_peri) (struct reb_simulation* const r, const unsigned int j) = ri_trace->S_peri ? ri_trace->S_peri : reb_integrator_trace_switch_peri_default; int new_close_encounter = 0; // New CEs // Clear encounter maps diff --git a/src/rebound.h b/src/rebound.h index fdf4585b6..f3431e2a0 100644 --- a/src/rebound.h +++ b/src/rebound.h @@ -849,7 +849,7 @@ DLLEXPORT double reb_integrator_mercurius_L_C5(const struct reb_simulation* cons // Built in trace switching functions -DLLEXPORT int reb_integrator_trace_switch_peri_pham(struct reb_simulation* const r, const unsigned int j); +DLLEXPORT int reb_integrator_trace_switch_peri_default(struct reb_simulation* const r, const unsigned int j); DLLEXPORT int reb_integrator_trace_switch_peri_fdot(struct reb_simulation* const r, const unsigned int j); DLLEXPORT int reb_integrator_trace_switch_peri_distance(struct reb_simulation* const r, const unsigned int j); DLLEXPORT int reb_integrator_trace_switch_peri_none(struct reb_simulation* const r, const unsigned int j);