From 42deafeb9490f9d98d1d7d9c292caa62bc7628b1 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Thu, 18 Jul 2024 16:11:11 -0400 Subject: [PATCH 01/16] encounter predict switching condition --- rebound/integrators/trace.py | 2 + src/integrator_trace.c | 144 +++++++++++++++++++++++++++++++++++ src/rebound.h | 5 +- 3 files changed, 150 insertions(+), 1 deletion(-) diff --git a/rebound/integrators/trace.py b/rebound/integrators/trace.py index 6f979be7d..a848db130 100644 --- a/rebound/integrators/trace.py +++ b/rebound/integrators/trace.py @@ -60,6 +60,8 @@ def __repr__(self): ("_particles_backup", ctypes.POINTER(Particle)), ("_particles_backup_kepler", ctypes.POINTER(Particle)), ("_particles_backup_additional_forces", ctypes.POINTER(Particle)), + ("_particles_pre", ctypes.POINTER(Particle)), + ("_particles_post", ctypes.POINTER(Particle)), ("_encounter_map", ctypes.POINTER(ctypes.c_int)), ("_com_pos", Vec3dBasic), ("_com_vel", Vec3dBasic), diff --git a/src/integrator_trace.c b/src/integrator_trace.c index 2194db43e..b6c5ad6e7 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -75,6 +75,96 @@ int reb_integrator_trace_switch_default(struct reb_simulation* const r, const un return d2*d2*d2 < dcritmax6; } +int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r, const unsigned int i, const unsigned int j){ + const struct reb_integrator_trace* const ri_trace = &(r->ri_trace); + const double dt = r->dt; + + + const double dx0 = ri_trace->particles_pre[i].x - ri_trace->particles_pre[j].x; + const double dy0 = ri_trace->particles_pre[i].y - ri_trace->particles_pre[j].y; + const double dz0 = ri_trace->particles_pre[i].z - ri_trace->particles_pre[j].z; + const double r0 = (dx0*dx0 + dy0*dy0 + dz0*dz0); + const double dvx0 = ri_trace->particles_pre[i].vx - ri_trace->particles_pre[j].vx; + const double dvy0 = ri_trace->particles_pre[i].vy - ri_trace->particles_pre[j].vy; + const double dvz0 = ri_trace->particles_pre[i].vz - ri_trace->particles_pre[j].vz; + const double dxn = ri_trace->particles_post[i].x - ri_trace->particles_post[j].x; + const double dyn = ri_trace->particles_post[i].y - ri_trace->particles_post[j].y; + const double dzn = ri_trace->particles_post[i].z - ri_trace->particles_post[j].z; + const double dvxn = ri_trace->particles_post[i].vx - ri_trace->particles_post[j].vx; + const double dvyn = ri_trace->particles_post[i].vy - ri_trace->particles_post[j].vy; + const double dvzn = ri_trace->particles_post[i].vz - ri_trace->particles_post[j].vz; + const double rn = (dxn*dxn + dyn*dyn + dzn*dzn); + + //const double dxp = r->particles[i].x - r->particles[j].x; + //const double dyp = r->particles[i].y - r->particles[j].y; + //const double dzp = r->particles[i].z - r->particles[j].z; + //const double rp = (dxp*dxp + dyp*dyp + dzp*dzp); + + + const double drndt = (dxn*dvxn+dyn*dvyn+dzn*dvzn)*2.; + const double drodt = (dx0*dvx0+dy0*dvy0+dz0*dvz0)*2.; + + const double a = 6.*(r0-rn)+3.*dt*(drodt+drndt); + const double b = 6.*(rn-r0)-2.*dt*(2.*drodt+drndt); + const double c = dt*drodt; + + double rmin = MIN(rn,r0); + + const double s = b*b-4.*a*c; + const double sr = sqrt(MAX(0.,s)); + const double tmin1 = (-b + sr)/(2.*a); + const double tmin2 = (-b - sr)/(2.*a); + if (tmin1>0. && tmin1<1.){ + const double rmin1 = (1.-tmin1)*(1.-tmin1)*(1.+2.*tmin1)*r0 + + tmin1*tmin1*(3.-2.*tmin1)*rn + + tmin1*(1.-tmin1)*(1.-tmin1)*dt*drodt + - tmin1*tmin1*(1.-tmin1)*dt*drndt; + rmin = MIN(MAX(rmin1,0.),rmin); + } + if (tmin2>0. && tmin2<1.){ + const double rmin2 = (1.-tmin2)*(1.-tmin2)*(1.+2.*tmin2)*r0 + + tmin2*tmin2*(3.-2.*tmin2)*rn + + tmin2*(1.-tmin2)*(1.-tmin2)*dt*drodt + - tmin2*tmin2*(1.-tmin2)*dt*drndt; + rmin = MIN(MAX(rmin2,0.),rmin); + } + + double dcriti6 = 0.0; + double dcritj6 = 0.0; + + const double m0 = r->particles[0].m; + + if (r->particles[i].m != 0){ + const double dxi = r->particles[i].x; // in dh + const double dyi = r->particles[i].y; + const double dzi = r->particles[i].z; + const double di2 = dxi*dxi + dyi*dyi + dzi*dzi; + const double mr = r->particles[i].m/(3.*m0); + dcriti6 = di2*di2*di2*mr*mr; + } + + if (r->particles[j].m != 0){ + const double dxj = r->particles[j].x; // in dh + const double dyj = r->particles[j].y; + const double dzj = r->particles[j].z; + const double dj2 = dxj*dxj + dyj*dyj + dzj*dzj; + const double mr = r->particles[j].m/(3.*m0); + dcritj6 = dj2*dj2*dj2*mr*mr; + } + + //const double dx = r->particles[i].x - r->particles[j].x; + //const double dy = r->particles[i].y - r->particles[j].y; + //const double dz = r->particles[i].z - r->particles[j].z; + //const double d2 = dx*dx + dy*dy + dz*dz; + //rmin = dx * dx + dy * dy + dz * dz; + + double r_crit_hill2 = ri_trace->r_crit_hill*ri_trace->r_crit_hill; + double dcritmax6 = r_crit_hill2 * r_crit_hill2 * r_crit_hill2 * MAX(dcriti6,dcritj6); + + return rmin*rmin*rmin < dcritmax6; + +} + int reb_integrator_trace_switch_peri_distance(struct reb_simulation* const r, const unsigned int j){ const struct reb_integrator_trace* const ri_trace = &(r->ri_trace); const double peri = ri_trace->peri_crit_distance; @@ -164,6 +254,53 @@ int reb_integrator_trace_switch_peri_default(struct reb_simulation* const r, con } } +int reb_integrator_trace_switch_peri_ep_accels(struct reb_simulation* const r, const unsigned int j){ + 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. + double h2 = r->dt/2.; + + double x = r->particles[j].x; + double y = r->particles[j].y; + double z = r->particles[j].z; + double d2 = x*x + y*y + z*z; + double d = sqrt(d2); + + // first derivative + double dx = r->particles[j].vx; + double dy = r->particles[j].vy; + double dz = r->particles[j].vz; + + // second derivative + double prefact2 = -GM/(d2*d); + double ddx = prefact2*x; + double ddy = prefact2*y; + double ddz = prefact2*z; + + // TLUEPP + + struct reb_particle pi_pre = {}; + pi_pre.x = r->particles[j].x - h2 * r->particles[j].vx; + pi_pre.y = r->particles[j].y - h2 * r->particles[j].vy; + pi_pre.z = r->particles[j].z - h2 * r->particles[j].vz; + pi_pre.vx = r->particles[j].vx - h2 * ddx; + pi_pre.vy = r->particles[j].vy - h2 * ddy; + pi_pre.vz = r->particles[j].vz - h2 * ddz; + + struct reb_particle pi_post = {}; + pi_post.x = r->particles[j].x + h2 * r->particles[j].vx; + pi_post.y = r->particles[j].y + h2 * r->particles[j].vy; + pi_post.z = r->particles[j].z + h2 * r->particles[j].vz; + pi_post.vx = r->particles[j].vx + h2 * ddx; + pi_post.vy = r->particles[j].vy + h2 * ddy; + pi_post.vz = r->particles[j].vz + h2 * ddz; + + ri_trace->particles_pre[j] = pi_pre; + ri_trace->particles_post[j] = pi_post; + + return 0; + +} + int reb_integrator_trace_switch_peri_none(struct reb_simulation* const r, const unsigned int j){ // No pericenter flags return 0; @@ -519,6 +656,8 @@ void reb_integrator_trace_part1(struct reb_simulation* r){ // These arrays are only used within one timestep. // Can be recreated without loosing bit-wise reproducibility. ri_trace->particles_backup = realloc(ri_trace->particles_backup,sizeof(struct reb_particle)*N); + ri_trace->particles_pre = realloc(ri_trace->particles_pre,sizeof(struct reb_particle)*N); + ri_trace->particles_post = realloc(ri_trace->particles_post,sizeof(struct reb_particle)*N); ri_trace->current_Ks = realloc(ri_trace->current_Ks,sizeof(int)*N*N); ri_trace->encounter_map = realloc(ri_trace->encounter_map,sizeof(int)*N); ri_trace->particles_backup_kepler = realloc(ri_trace->particles_backup_kepler,sizeof(struct reb_particle)*N); @@ -841,6 +980,11 @@ void reb_integrator_trace_reset(struct reb_simulation* r){ r->ri_trace.particles_backup_kepler = NULL; free(r->ri_trace.particles_backup_additional_forces); r->ri_trace.particles_backup_additional_forces = NULL; + + free(r->ri_trace.particles_pre); + r->ri_trace.particles_pre = NULL; + free(r->ri_trace.particles_post); + r->ri_trace.particles_post = NULL; free(r->ri_trace.encounter_map); diff --git a/src/rebound.h b/src/rebound.h index e1f239697..284e29f98 100644 --- a/src/rebound.h +++ b/src/rebound.h @@ -291,6 +291,8 @@ struct reb_integrator_trace { struct reb_particle* REB_RESTRICT particles_backup; // Contains coordinates before the entire step struct reb_particle* REB_RESTRICT particles_backup_kepler; // Contains coordinates before kepler step struct reb_particle* REB_RESTRICT particles_backup_additional_forces; // For additional forces + struct reb_particle* REB_RESTRICT particles_pre; // Pre-Timestep coordinates + struct reb_particle* REB_RESTRICT particles_post; // Post-Timestep coordinates int* encounter_map; // Map to represent which particles are integrated with BS struct reb_vec3d com_pos; // Used to keep track of the centre of mass during the timestep @@ -849,11 +851,12 @@ DLLEXPORT double reb_integrator_mercurius_L_C5(const struct reb_simulation* cons // Built in trace switching functions -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); +DLLEXPORT int reb_integrator_trace_switch_peri_ep_accels(struct reb_simulation* const r, const unsigned int j); // This is a placeholder DLLEXPORT int reb_integrator_trace_switch_default(struct reb_simulation* const r, const unsigned int i, const unsigned int j); +DLLEXPORT int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r, const unsigned int i, const unsigned int j); // Built in collision resolve functions From 61398c001b7448f875fcb37da6612af4c81c520f Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Thu, 18 Jul 2024 16:26:03 -0400 Subject: [PATCH 02/16] typo --- src/integrator_trace.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index b6c5ad6e7..a42f74a35 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -152,11 +152,11 @@ int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r dcritj6 = dj2*dj2*dj2*mr*mr; } - //const double dx = r->particles[i].x - r->particles[j].x; - //const double dy = r->particles[i].y - r->particles[j].y; - //const double dz = r->particles[i].z - r->particles[j].z; + const double dx = r->particles[i].x - r->particles[j].x; + const double dy = r->particles[i].y - r->particles[j].y; + const double dz = r->particles[i].z - r->particles[j].z; //const double d2 = dx*dx + dy*dy + dz*dz; - //rmin = dx * dx + dy * dy + dz * dz; + rmin = dx * dx + dy * dy + dz * dz; double r_crit_hill2 = ri_trace->r_crit_hill*ri_trace->r_crit_hill; double dcritmax6 = r_crit_hill2 * r_crit_hill2 * r_crit_hill2 * MAX(dcriti6,dcritj6); From fbf4c5ca2276f913b7030c0ad6a2f4bb89c9a12a Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Thu, 18 Jul 2024 19:48:00 -0400 Subject: [PATCH 03/16] fix empty particle syntax --- src/integrator_trace.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index a42f74a35..93ccf4b62 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -278,7 +278,7 @@ int reb_integrator_trace_switch_peri_ep_accels(struct reb_simulation* const r, c // TLUEPP - struct reb_particle pi_pre = {}; + struct reb_particle pi_pre = {0}; pi_pre.x = r->particles[j].x - h2 * r->particles[j].vx; pi_pre.y = r->particles[j].y - h2 * r->particles[j].vy; pi_pre.z = r->particles[j].z - h2 * r->particles[j].vz; @@ -286,7 +286,7 @@ int reb_integrator_trace_switch_peri_ep_accels(struct reb_simulation* const r, c pi_pre.vy = r->particles[j].vy - h2 * ddy; pi_pre.vz = r->particles[j].vz - h2 * ddz; - struct reb_particle pi_post = {}; + struct reb_particle pi_post = {0}; pi_post.x = r->particles[j].x + h2 * r->particles[j].vx; pi_post.y = r->particles[j].y + h2 * r->particles[j].vy; pi_post.z = r->particles[j].z + h2 * r->particles[j].vz; From 3ee18ad79c92475d905e1836c1f9069be127d77e Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Thu, 18 Jul 2024 20:35:19 -0400 Subject: [PATCH 04/16] added acceleration term --- src/integrator_trace.c | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index 93ccf4b62..0c35055c4 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -152,11 +152,11 @@ int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r dcritj6 = dj2*dj2*dj2*mr*mr; } - const double dx = r->particles[i].x - r->particles[j].x; - const double dy = r->particles[i].y - r->particles[j].y; - const double dz = r->particles[i].z - r->particles[j].z; +// const double dx = r->particles[i].x - r->particles[j].x; +// const double dy = r->particles[i].y - r->particles[j].y; +// const double dz = r->particles[i].z - r->particles[j].z; //const double d2 = dx*dx + dy*dy + dz*dz; - rmin = dx * dx + dy * dy + dz * dz; +// rmin = dx * dx + dy * dy + dz * dz; double r_crit_hill2 = ri_trace->r_crit_hill*ri_trace->r_crit_hill; double dcritmax6 = r_crit_hill2 * r_crit_hill2 * r_crit_hill2 * MAX(dcriti6,dcritj6); @@ -279,17 +279,17 @@ int reb_integrator_trace_switch_peri_ep_accels(struct reb_simulation* const r, c // TLUEPP struct reb_particle pi_pre = {0}; - pi_pre.x = r->particles[j].x - h2 * r->particles[j].vx; - pi_pre.y = r->particles[j].y - h2 * r->particles[j].vy; - pi_pre.z = r->particles[j].z - h2 * r->particles[j].vz; + pi_pre.x = r->particles[j].x - h2 * r->particles[j].vx - 0.5 * h2 * h2 * ddx; + pi_pre.y = r->particles[j].y - h2 * r->particles[j].vy - 0.5 * h2 * h2 * ddy; + pi_pre.z = r->particles[j].z - h2 * r->particles[j].vz - 0.5 * h2 * h2 * ddz; pi_pre.vx = r->particles[j].vx - h2 * ddx; pi_pre.vy = r->particles[j].vy - h2 * ddy; pi_pre.vz = r->particles[j].vz - h2 * ddz; struct reb_particle pi_post = {0}; - pi_post.x = r->particles[j].x + h2 * r->particles[j].vx; - pi_post.y = r->particles[j].y + h2 * r->particles[j].vy; - pi_post.z = r->particles[j].z + h2 * r->particles[j].vz; + pi_post.x = r->particles[j].x + h2 * r->particles[j].vx + 0.5 * h2 * h2 * ddx; + pi_post.y = r->particles[j].y + h2 * r->particles[j].vy + 0.5 * h2 * h2 * ddy; + pi_post.z = r->particles[j].z + h2 * r->particles[j].vz + 0.5 * h2 * h2 * ddz; pi_post.vx = r->particles[j].vx + h2 * ddx; pi_post.vy = r->particles[j].vy + h2 * ddy; pi_post.vz = r->particles[j].vz + h2 * ddz; From 6b4c9192998cc250e2319a32e6792d0124bdd156 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Sat, 20 Jul 2024 14:25:35 -0400 Subject: [PATCH 05/16] added line encounter predict --- src/integrator_trace.c | 146 ++++++++++++++++++++++++++++------------- src/rebound.h | 1 + 2 files changed, 102 insertions(+), 45 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index 0c35055c4..79b8f330f 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -80,54 +80,46 @@ int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r const double dt = r->dt; + // pre const double dx0 = ri_trace->particles_pre[i].x - ri_trace->particles_pre[j].x; const double dy0 = ri_trace->particles_pre[i].y - ri_trace->particles_pre[j].y; const double dz0 = ri_trace->particles_pre[i].z - ri_trace->particles_pre[j].z; const double r0 = (dx0*dx0 + dy0*dy0 + dz0*dz0); - const double dvx0 = ri_trace->particles_pre[i].vx - ri_trace->particles_pre[j].vx; - const double dvy0 = ri_trace->particles_pre[i].vy - ri_trace->particles_pre[j].vy; - const double dvz0 = ri_trace->particles_pre[i].vz - ri_trace->particles_pre[j].vz; + //const double dvx0 = ri_trace->particles_pre[i].vx - ri_trace->particles_pre[j].vx; + //const double dvy0 = ri_trace->particles_pre[i].vy - ri_trace->particles_pre[j].vy; + //const double dvz0 = ri_trace->particles_pre[i].vz - ri_trace->particles_pre[j].vz; + + // post const double dxn = ri_trace->particles_post[i].x - ri_trace->particles_post[j].x; const double dyn = ri_trace->particles_post[i].y - ri_trace->particles_post[j].y; const double dzn = ri_trace->particles_post[i].z - ri_trace->particles_post[j].z; - const double dvxn = ri_trace->particles_post[i].vx - ri_trace->particles_post[j].vx; - const double dvyn = ri_trace->particles_post[i].vy - ri_trace->particles_post[j].vy; - const double dvzn = ri_trace->particles_post[i].vz - ri_trace->particles_post[j].vz; + //const double dvxn = ri_trace->particles_post[i].vx - ri_trace->particles_post[j].vx; + //const double dvyn = ri_trace->particles_post[i].vy - ri_trace->particles_post[j].vy; + //const double dvzn = ri_trace->particles_post[i].vz - ri_trace->particles_post[j].vz; const double rn = (dxn*dxn + dyn*dyn + dzn*dzn); - //const double dxp = r->particles[i].x - r->particles[j].x; - //const double dyp = r->particles[i].y - r->particles[j].y; - //const double dzp = r->particles[i].z - r->particles[j].z; - //const double rp = (dxp*dxp + dyp*dyp + dzp*dzp); - - - const double drndt = (dxn*dvxn+dyn*dvyn+dzn*dvzn)*2.; - const double drodt = (dx0*dvx0+dy0*dvy0+dz0*dvz0)*2.; - - const double a = 6.*(r0-rn)+3.*dt*(drodt+drndt); - const double b = 6.*(rn-r0)-2.*dt*(2.*drodt+drndt); - const double c = dt*drodt; - - double rmin = MIN(rn,r0); - - const double s = b*b-4.*a*c; - const double sr = sqrt(MAX(0.,s)); - const double tmin1 = (-b + sr)/(2.*a); - const double tmin2 = (-b - sr)/(2.*a); - if (tmin1>0. && tmin1<1.){ - const double rmin1 = (1.-tmin1)*(1.-tmin1)*(1.+2.*tmin1)*r0 - + tmin1*tmin1*(3.-2.*tmin1)*rn - + tmin1*(1.-tmin1)*(1.-tmin1)*dt*drodt - - tmin1*tmin1*(1.-tmin1)*dt*drndt; - rmin = MIN(MAX(rmin1,0.),rmin); - } - if (tmin2>0. && tmin2<1.){ - const double rmin2 = (1.-tmin2)*(1.-tmin2)*(1.+2.*tmin2)*r0 - + tmin2*tmin2*(3.-2.*tmin2)*rn - + tmin2*(1.-tmin2)*(1.-tmin2)*dt*drodt - - tmin2*tmin2*(1.-tmin2)*dt*drndt; - rmin = MIN(MAX(rmin2,0.),rmin); - } + // current + const double dxi = r->particles[i].x; // in dh + const double dyi = r->particles[i].y; + const double dzi = r->particles[i].z; + const double dxj = r->particles[j].x; // in dh + const double dyj = r->particles[j].y; + const double dzj = r->particles[j].z; + const double dxp = dxi - dxj; + const double dyp = dyi - dyj; + const double dzp = dzi - dzj; + const double rp = (dxp*dxp + dyp*dyp + dzp*dzp); + + //const double drndt = (dxn*dvxn+dyn*dvyn+dzn*dvzn)*2.; + //const double drodt = (dx0*dvx0+dy0*dvy0+dz0*dvz0)*2.; + + const double a = -2.*(2.*rp - rn - r0)/dt; + const double b = (rn-rp)/dt; + const double c = rp; + + const double v = -b/(2*a); + const double vmin = a*v*v + b*v + c; + double rmin = MIN(MIN(MIN(rn,r0),rp),vmin*vmin); double dcriti6 = 0.0; double dcritj6 = 0.0; @@ -135,18 +127,12 @@ int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r const double m0 = r->particles[0].m; if (r->particles[i].m != 0){ - const double dxi = r->particles[i].x; // in dh - const double dyi = r->particles[i].y; - const double dzi = r->particles[i].z; const double di2 = dxi*dxi + dyi*dyi + dzi*dzi; const double mr = r->particles[i].m/(3.*m0); dcriti6 = di2*di2*di2*mr*mr; } if (r->particles[j].m != 0){ - const double dxj = r->particles[j].x; // in dh - const double dyj = r->particles[j].y; - const double dzj = r->particles[j].z; const double dj2 = dxj*dxj + dyj*dyj + dzj*dzj; const double mr = r->particles[j].m/(3.*m0); dcritj6 = dj2*dj2*dj2*mr*mr; @@ -165,6 +151,76 @@ int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r } +int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, const unsigned int i, const unsigned int j){ + struct reb_integrator_trace* const ri_trace = &(r->ri_trace); + const double h2 = r->dt/2.; + + const double dxi = r->particles[i].x; + const double dyi = r->particles[i].y; + const double dzi = r->particles[i].z; + + const double dxj = r->particles[j].x; + const double dyj = r->particles[j].y; + const double dzj = r->particles[j].z; + + const double dx = dxi - dxj; + const double dy = dyi - dyj; + const double dz = dzi - dzj; + const double rp = dx*dx + dy*dy + dz*dz; + + const double dvx = r->particles[i].vx - r->particles[j].vx; + const double dvy = r->particles[i].vy - r->particles[j].vy; + const double dvz = r->particles[i].vz - r->particles[j].vz; + + // pre + const double dx_pre = dx - h2 * dvx; + const double dy_pre = dy - h2 * dvy; + const double dz_pre = dz - h2 * dvz; + const double rpre = dx_pre*dx_pre + dy_pre*dy_pre + dz_pre*dz_pre; + + // post + const double dx_post = dx + h2 * dvx; + const double dy_post = dy + h2 * dvy; + const double dz_post = dz + h2 * dvz; + const double rpost = dx_post*dx_post + dy_post*dy_post + dz_post*dz_post; + + double rmin1 = MIN(rpre,rpost); + double rmin = MIN(rmin1,rp); + + const double a = dvx*dvx + dvy*dvy + dvz*dvz; + const double b = 2.*(dx*dvx+dy*dvy+dz*dvz); + const double c = dx*dx + dy*dy + dz*dz; + + const double tmin = -b/(2*a); + if (tmin > -1.*h2 && tmin < h2){ + double rmin2 = a*tmin*tmin + b*tmin + c; + rmin = MIN(rmin2, rmin); + } + + + double dcriti6 = 0.0; + double dcritj6 = 0.0; + + const double m0 = r->particles[0].m; + + if (r->particles[i].m != 0){ + const double di2 = dxi*dxi + dyi*dyi + dzi*dzi; + const double mr = r->particles[i].m/(3.*m0); + dcriti6 = di2*di2*di2*mr*mr; + } + + if (r->particles[j].m != 0){ + const double dj2 = dxj*dxj + dyj*dyj + dzj*dzj; + const double mr = r->particles[j].m/(3.*m0); + dcritj6 = dj2*dj2*dj2*mr*mr; + } + + double r_crit_hill2 = ri_trace->r_crit_hill*ri_trace->r_crit_hill; + double dcritmax6 = r_crit_hill2 * r_crit_hill2 * r_crit_hill2 * MAX(dcriti6,dcritj6); + + return rmin*rmin*rmin < dcritmax6; +} + int reb_integrator_trace_switch_peri_distance(struct reb_simulation* const r, const unsigned int j){ const struct reb_integrator_trace* const ri_trace = &(r->ri_trace); const double peri = ri_trace->peri_crit_distance; diff --git a/src/rebound.h b/src/rebound.h index 284e29f98..d2e86cfea 100644 --- a/src/rebound.h +++ b/src/rebound.h @@ -857,6 +857,7 @@ DLLEXPORT int reb_integrator_trace_switch_peri_none(struct reb_simulation* const DLLEXPORT int reb_integrator_trace_switch_peri_ep_accels(struct reb_simulation* const r, const unsigned int j); // This is a placeholder DLLEXPORT int reb_integrator_trace_switch_default(struct reb_simulation* const r, const unsigned int i, const unsigned int j); DLLEXPORT int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r, const unsigned int i, const unsigned int j); +DLLEXPORT int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, const unsigned int i, const unsigned int j); // Built in collision resolve functions From 5da110b520605537a9b11f6ef1af44b6092f6669 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Sat, 20 Jul 2024 17:30:14 -0400 Subject: [PATCH 06/16] code cleanup, cubic interpolation and unit test for encounter prediction --- rebound/integrators/trace.py | 7 +- rebound/tests/test_trace_basic.py | 14 +++ src/integrator_trace.c | 137 ++++++++++++------------------ src/rebound.h | 3 - 4 files changed, 72 insertions(+), 89 deletions(-) diff --git a/rebound/integrators/trace.py b/rebound/integrators/trace.py index a848db130..104ddbfe2 100644 --- a/rebound/integrators/trace.py +++ b/rebound/integrators/trace.py @@ -60,8 +60,6 @@ def __repr__(self): ("_particles_backup", ctypes.POINTER(Particle)), ("_particles_backup_kepler", ctypes.POINTER(Particle)), ("_particles_backup_additional_forces", ctypes.POINTER(Particle)), - ("_particles_pre", ctypes.POINTER(Particle)), - ("_particles_post", ctypes.POINTER(Particle)), ("_encounter_map", ctypes.POINTER(ctypes.c_int)), ("_com_pos", Vec3dBasic), ("_com_vel", Vec3dBasic), @@ -69,7 +67,6 @@ def __repr__(self): ("_current_C", ctypes.c_uint), ("_force_accept", ctypes.c_uint), ] - # To be honest I'm not sure what these do: do we need this? - Tiger @property def S(self): raise AttributeError("You can only set C function pointers from python.") @@ -77,6 +74,10 @@ def S(self): def S(self, func): if func == "default": self._S = cast(clibrebound.reb_integrator_trace_switch_default,TRACEKF) + elif func == "line": + self._S = cast(clibrebound.reb_integrator_trace_switch_encounter_line, TRACEKF) + elif func == "cubic": + self._S = cast(clibrebound.reb_integrator_trace_switch_encounter_predict, TRACEKF) else: self._Sfp = TRACEKF(func) self._S = self._Sfp diff --git a/rebound/tests/test_trace_basic.py b/rebound/tests/test_trace_basic.py index 604a1366c..ec8991a5a 100644 --- a/rebound/tests/test_trace_basic.py +++ b/rebound/tests/test_trace_basic.py @@ -36,6 +36,20 @@ def test_trace_encounter_condition(self): for j in range(i+1,sim.N): self.assertEqual(sim.ri_trace._current_Ks[i*sim.N+j],0) + def test_trace_encounter_prediction(self): + sim = rebound.Simulation() + sim.add(m=1) + sim.add(m=9.55e-4,x=5.2) + sim.add(x=5.3,y=0.36,vy=-7.2) # Normal switching condition misses this + sim.integrator = "TRACE" + sim.dt = 0.01 + sim.ri_trace.r_crit_hill = 1 + sim.S = "line" + sim.step() + + self.assertEqual(sim.ri_trace._encounter_N,3) + + if __name__ == "__main__": unittest.main() diff --git a/src/integrator_trace.c b/src/integrator_trace.c index 79b8f330f..4badfcb6d 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -78,25 +78,9 @@ int reb_integrator_trace_switch_default(struct reb_simulation* const r, const un int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r, const unsigned int i, const unsigned int j){ const struct reb_integrator_trace* const ri_trace = &(r->ri_trace); const double dt = r->dt; - - - // pre - const double dx0 = ri_trace->particles_pre[i].x - ri_trace->particles_pre[j].x; - const double dy0 = ri_trace->particles_pre[i].y - ri_trace->particles_pre[j].y; - const double dz0 = ri_trace->particles_pre[i].z - ri_trace->particles_pre[j].z; - const double r0 = (dx0*dx0 + dy0*dy0 + dz0*dz0); - //const double dvx0 = ri_trace->particles_pre[i].vx - ri_trace->particles_pre[j].vx; - //const double dvy0 = ri_trace->particles_pre[i].vy - ri_trace->particles_pre[j].vy; - //const double dvz0 = ri_trace->particles_pre[i].vz - ri_trace->particles_pre[j].vz; + const double h2 = dt/2.; - // post - const double dxn = ri_trace->particles_post[i].x - ri_trace->particles_post[j].x; - const double dyn = ri_trace->particles_post[i].y - ri_trace->particles_post[j].y; - const double dzn = ri_trace->particles_post[i].z - ri_trace->particles_post[j].z; - //const double dvxn = ri_trace->particles_post[i].vx - ri_trace->particles_post[j].vx; - //const double dvyn = ri_trace->particles_post[i].vy - ri_trace->particles_post[j].vy; - //const double dvzn = ri_trace->particles_post[i].vz - ri_trace->particles_post[j].vz; - const double rn = (dxn*dxn + dyn*dyn + dzn*dzn); + double GM = r->G*r->particles[0].m; // Think this is still the correct mass, for these encounters solar mass should be more relevant // current const double dxi = r->particles[i].x; // in dh @@ -109,17 +93,64 @@ int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r const double dyp = dyi - dyj; const double dzp = dzi - dzj; const double rp = (dxp*dxp + dyp*dyp + dzp*dzp); + + const double di2 = dxi*dxi + dyi*dyi + dzi*dzi; + const double di = sqrt(di2); + const double dj2 = dxj*dxj + dyj*dyj + dzj*dzj; + const double dj = sqrt(dj2); + + // second derivative + double prefacti2 = -GM/(di2*di); + double ddxi = prefacti2*dxi; + double ddyi = prefacti2*dyi; + double ddzi = prefacti2*dzi; + + double prefactj2 = -GM/(dj2*dj); + double ddxj = prefactj2*dxj; + double ddyj = prefactj2*dyj; + double ddzj = prefactj2*dzj; + + // TLUEPP + + const double pi_pre_x = r->particles[i].x - h2 * r->particles[i].vx + 0.5 * h2 * h2 * ddxi; + const double pi_pre_y = r->particles[i].y - h2 * r->particles[i].vy + 0.5 * h2 * h2 * ddyi; + const double pi_pre_z = r->particles[i].z - h2 * r->particles[i].vz + 0.5 * h2 * h2 * ddzi; + + const double pi_post_x = r->particles[i].x + h2 * r->particles[i].vx + 0.5 * h2 * h2 * ddxi; + const double pi_post_y = r->particles[i].y + h2 * r->particles[i].vy + 0.5 * h2 * h2 * ddyi; + const double pi_post_z = r->particles[i].z + h2 * r->particles[i].vz + 0.5 * h2 * h2 * ddzi; - //const double drndt = (dxn*dvxn+dyn*dvyn+dzn*dvzn)*2.; - //const double drodt = (dx0*dvx0+dy0*dvy0+dz0*dvz0)*2.; + const double pj_pre_x = r->particles[j].x - h2 * r->particles[j].vx + 0.5 * h2 * h2 * ddxj; + const double pj_pre_y = r->particles[j].y - h2 * r->particles[j].vy + 0.5 * h2 * h2 * ddyj; + const double pj_pre_z = r->particles[j].z - h2 * r->particles[j].vz + 0.5 * h2 * h2 * ddzj; + + const double pj_post_x = r->particles[j].x + h2 * r->particles[j].vx + 0.5 * h2 * h2 * ddxj; + const double pj_post_y = r->particles[j].y + h2 * r->particles[j].vy + 0.5 * h2 * h2 * ddyj; + const double pj_post_z = r->particles[j].z + h2 * r->particles[j].vz + 0.5 * h2 * h2 * ddzj; + + // pre + const double dx0 = pi_pre_x - pj_pre_x; + const double dy0 = pi_pre_y - pj_pre_y; + const double dz0 = pi_pre_z - pj_pre_z; + const double r0 = (dx0*dx0 + dy0*dy0 + dz0*dz0); + + // post + const double dxn = pi_post_x - pj_post_x; + const double dyn = pi_post_y - pj_post_y; + const double dzn = pi_post_z - pj_post_z; + const double rn = (dxn*dxn + dyn*dyn + dzn*dzn); const double a = -2.*(2.*rp - rn - r0)/dt; const double b = (rn-rp)/dt; const double c = rp; + double rmin = MIN(MIN(rn,r0),rp); const double v = -b/(2*a); - const double vmin = a*v*v + b*v + c; - double rmin = MIN(MIN(MIN(rn,r0),rp),vmin*vmin); + if (v > -1.*h2 && v < h2){ + const double vmin = a*v*v + b*v + c; + rmin = MIN(rmin, vmin*vmin); + } + double dcriti6 = 0.0; double dcritj6 = 0.0; @@ -138,11 +169,6 @@ int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r dcritj6 = dj2*dj2*dj2*mr*mr; } -// const double dx = r->particles[i].x - r->particles[j].x; -// const double dy = r->particles[i].y - r->particles[j].y; -// const double dz = r->particles[i].z - r->particles[j].z; - //const double d2 = dx*dx + dy*dy + dz*dz; -// rmin = dx * dx + dy * dy + dz * dz; double r_crit_hill2 = ri_trace->r_crit_hill*ri_trace->r_crit_hill; double dcritmax6 = r_crit_hill2 * r_crit_hill2 * r_crit_hill2 * MAX(dcriti6,dcritj6); @@ -310,53 +336,6 @@ int reb_integrator_trace_switch_peri_default(struct reb_simulation* const r, con } } -int reb_integrator_trace_switch_peri_ep_accels(struct reb_simulation* const r, const unsigned int j){ - 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. - double h2 = r->dt/2.; - - double x = r->particles[j].x; - double y = r->particles[j].y; - double z = r->particles[j].z; - double d2 = x*x + y*y + z*z; - double d = sqrt(d2); - - // first derivative - double dx = r->particles[j].vx; - double dy = r->particles[j].vy; - double dz = r->particles[j].vz; - - // second derivative - double prefact2 = -GM/(d2*d); - double ddx = prefact2*x; - double ddy = prefact2*y; - double ddz = prefact2*z; - - // TLUEPP - - struct reb_particle pi_pre = {0}; - pi_pre.x = r->particles[j].x - h2 * r->particles[j].vx - 0.5 * h2 * h2 * ddx; - pi_pre.y = r->particles[j].y - h2 * r->particles[j].vy - 0.5 * h2 * h2 * ddy; - pi_pre.z = r->particles[j].z - h2 * r->particles[j].vz - 0.5 * h2 * h2 * ddz; - pi_pre.vx = r->particles[j].vx - h2 * ddx; - pi_pre.vy = r->particles[j].vy - h2 * ddy; - pi_pre.vz = r->particles[j].vz - h2 * ddz; - - struct reb_particle pi_post = {0}; - pi_post.x = r->particles[j].x + h2 * r->particles[j].vx + 0.5 * h2 * h2 * ddx; - pi_post.y = r->particles[j].y + h2 * r->particles[j].vy + 0.5 * h2 * h2 * ddy; - pi_post.z = r->particles[j].z + h2 * r->particles[j].vz + 0.5 * h2 * h2 * ddz; - pi_post.vx = r->particles[j].vx + h2 * ddx; - pi_post.vy = r->particles[j].vy + h2 * ddy; - pi_post.vz = r->particles[j].vz + h2 * ddz; - - ri_trace->particles_pre[j] = pi_pre; - ri_trace->particles_post[j] = pi_post; - - return 0; - -} - int reb_integrator_trace_switch_peri_none(struct reb_simulation* const r, const unsigned int j){ // No pericenter flags return 0; @@ -712,11 +691,9 @@ void reb_integrator_trace_part1(struct reb_simulation* r){ // These arrays are only used within one timestep. // Can be recreated without loosing bit-wise reproducibility. ri_trace->particles_backup = realloc(ri_trace->particles_backup,sizeof(struct reb_particle)*N); - ri_trace->particles_pre = realloc(ri_trace->particles_pre,sizeof(struct reb_particle)*N); - ri_trace->particles_post = realloc(ri_trace->particles_post,sizeof(struct reb_particle)*N); + ri_trace->particles_backup_kepler = realloc(ri_trace->particles_backup_kepler,sizeof(struct reb_particle)*N); ri_trace->current_Ks = realloc(ri_trace->current_Ks,sizeof(int)*N*N); ri_trace->encounter_map = realloc(ri_trace->encounter_map,sizeof(int)*N); - ri_trace->particles_backup_kepler = realloc(ri_trace->particles_backup_kepler,sizeof(struct reb_particle)*N); ri_trace->N_allocated = N; } @@ -1036,12 +1013,6 @@ void reb_integrator_trace_reset(struct reb_simulation* r){ r->ri_trace.particles_backup_kepler = NULL; free(r->ri_trace.particles_backup_additional_forces); r->ri_trace.particles_backup_additional_forces = NULL; - - free(r->ri_trace.particles_pre); - r->ri_trace.particles_pre = NULL; - free(r->ri_trace.particles_post); - r->ri_trace.particles_post = NULL; - free(r->ri_trace.encounter_map); r->ri_trace.encounter_map = NULL; diff --git a/src/rebound.h b/src/rebound.h index d2e86cfea..2d802b8e8 100644 --- a/src/rebound.h +++ b/src/rebound.h @@ -291,8 +291,6 @@ struct reb_integrator_trace { struct reb_particle* REB_RESTRICT particles_backup; // Contains coordinates before the entire step struct reb_particle* REB_RESTRICT particles_backup_kepler; // Contains coordinates before kepler step struct reb_particle* REB_RESTRICT particles_backup_additional_forces; // For additional forces - struct reb_particle* REB_RESTRICT particles_pre; // Pre-Timestep coordinates - struct reb_particle* REB_RESTRICT particles_post; // Post-Timestep coordinates int* encounter_map; // Map to represent which particles are integrated with BS struct reb_vec3d com_pos; // Used to keep track of the centre of mass during the timestep @@ -854,7 +852,6 @@ DLLEXPORT double reb_integrator_mercurius_L_C5(const struct reb_simulation* cons 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); -DLLEXPORT int reb_integrator_trace_switch_peri_ep_accels(struct reb_simulation* const r, const unsigned int j); // This is a placeholder DLLEXPORT int reb_integrator_trace_switch_default(struct reb_simulation* const r, const unsigned int i, const unsigned int j); DLLEXPORT int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r, const unsigned int i, const unsigned int j); DLLEXPORT int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, const unsigned int i, const unsigned int j); From a326a22066ae333b99b21429668905631d657998 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Tue, 23 Jul 2024 13:53:56 -0400 Subject: [PATCH 07/16] encounter predicts approximates as quartic, may have bugs --- src/integrator_trace.c | 80 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 71 insertions(+), 9 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index 4badfcb6d..6bd3a0fb1 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -115,18 +115,30 @@ int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r const double pi_pre_x = r->particles[i].x - h2 * r->particles[i].vx + 0.5 * h2 * h2 * ddxi; const double pi_pre_y = r->particles[i].y - h2 * r->particles[i].vy + 0.5 * h2 * h2 * ddyi; const double pi_pre_z = r->particles[i].z - h2 * r->particles[i].vz + 0.5 * h2 * h2 * ddzi; + const double pi_pre_vx = r->particles[i].vx - h2 * ddxi; + const double pi_pre_vy = r->particles[i].vy - h2 * ddyi; + const double pi_pre_vz = r->particles[i].vz - h2 * ddzi; const double pi_post_x = r->particles[i].x + h2 * r->particles[i].vx + 0.5 * h2 * h2 * ddxi; const double pi_post_y = r->particles[i].y + h2 * r->particles[i].vy + 0.5 * h2 * h2 * ddyi; const double pi_post_z = r->particles[i].z + h2 * r->particles[i].vz + 0.5 * h2 * h2 * ddzi; + const double pi_post_vx = r->particles[i].vx + h2 * ddxi; + const double pi_post_vy = r->particles[i].vy + h2 * ddyi; + const double pi_post_vz = r->particles[i].vz + h2 * ddzi; const double pj_pre_x = r->particles[j].x - h2 * r->particles[j].vx + 0.5 * h2 * h2 * ddxj; const double pj_pre_y = r->particles[j].y - h2 * r->particles[j].vy + 0.5 * h2 * h2 * ddyj; const double pj_pre_z = r->particles[j].z - h2 * r->particles[j].vz + 0.5 * h2 * h2 * ddzj; + const double pj_pre_vx = r->particles[j].vx - h2 * ddxj; + const double pj_pre_vy = r->particles[j].vy - h2 * ddyj; + const double pj_pre_vz = r->particles[j].vz - h2 * ddzj; const double pj_post_x = r->particles[j].x + h2 * r->particles[j].vx + 0.5 * h2 * h2 * ddxj; const double pj_post_y = r->particles[j].y + h2 * r->particles[j].vy + 0.5 * h2 * h2 * ddyj; const double pj_post_z = r->particles[j].z + h2 * r->particles[j].vz + 0.5 * h2 * h2 * ddzj; + const double pj_post_vx = r->particles[j].vx + h2 * ddxj; + const double pj_post_vy = r->particles[j].vy + h2 * ddyj; + const double pj_post_vz = r->particles[j].vz + h2 * ddzj; // pre const double dx0 = pi_pre_x - pj_pre_x; @@ -140,17 +152,67 @@ int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r const double dzn = pi_post_z - pj_post_z; const double rn = (dxn*dxn + dyn*dyn + dzn*dzn); - const double a = -2.*(2.*rp - rn - r0)/dt; - const double b = (rn-rp)/dt; - const double c = rp; - double rmin = MIN(MIN(rn,r0),rp); - - const double v = -b/(2*a); - if (v > -1.*h2 && v < h2){ - const double vmin = a*v*v + b*v + c; - rmin = MIN(rmin, vmin*vmin); + double rmin = MIN(MIN(rn, rp), r0); + + // Quartic Polynomial + const double a = 2.*(8.*rp + dt*(rn - r0) - 4.*(rn - r0)) / (dt*dt*dt*dt); + const double b = (dt*(rn + r0) - 2.*(rn - r0)) / (dt*dt*dt); + const double c = (8.*(rn + r0) - 16.*rp - dt*(rn-r0)) / (2*dt*dt); + const double d = (6.*(rn - r0) - dt*(rn + r0)) / (4. * dt); + const double e = rp; + + // Derivative -- cubic + const double ca = 4.*a; + const double cb = 3.*b; + const double cc = 2.*c; + const double cd = d; + + // Depressed cubic + const double p = (3.*ca*cc - cb*cb)/(3*ca*ca); + const double q = (2.*cb*cb*cb - 9.*ca*cb*cc + 27.*ca*ca*cd)/(27.*ca*ca*ca); + double roots[3] = {}; + unsigned int nroots = 0; // real roots + if (fabs(p) < 1e-10){ // p = 0 -> t^3 = -q + roots[0] = pow(-1.*q, 1./3.); + nroots = 1; + } + else if (fabs(q) < 1e-10 && p < 0){ // q = 0 -> t^3 + pt = 0 + roots[0] = 0.; + roots[1] = sqrt(-1.*p); + roots[2] = -1.*sqrt(-1.*p); + nroots = 3; + } + else{ + const double disc = q*q/4. + p*p*p/27.; + if (fabs(disc) < 1e-10){ + roots[0] = -1.5*p*q; + roots[1] = 3.*q/p; + nroots = 2; + } + else if (disc > 0.){ + const double u = pow(-q/2. - sqrt(d), 1./3.); + roots[0] = u - p/(3.*u); + nroots = 1; + } + else{ + const double u = 2.*sqrt(-p/3.); + const double t = acos(3.*q/p/u)/3.; + const double k = 2.*M_PI/3.; + roots[0] = u*cos(t); + roots[1] = u*cos(t-k); + roots[2] = u*cos(t-2*k); + nroots = 3; + } } + // check if roots are in the range + for (unsigned int i = 0; i < nroots; i++){ + double root = roots[i] - b/(3.*a); + if (root < h2 && root > -1.*h2){ + double rval = a*root*root*root*root + b*root*root*root + c*root*root + d*root + e; + rmin = MIN(rmin, rval); + } + } double dcriti6 = 0.0; double dcritj6 = 0.0; From b72a04a472ebcd1f805f7ea70b88ce54c6560af6 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Sat, 27 Jul 2024 10:24:56 -0400 Subject: [PATCH 08/16] reversibility fix and added pericenter debugging --- rebound/integrators/trace.py | 2 -- src/rebound.h | 3 ++- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/rebound/integrators/trace.py b/rebound/integrators/trace.py index 104ddbfe2..a78acb840 100644 --- a/rebound/integrators/trace.py +++ b/rebound/integrators/trace.py @@ -76,8 +76,6 @@ def S(self, func): self._S = cast(clibrebound.reb_integrator_trace_switch_default,TRACEKF) elif func == "line": self._S = cast(clibrebound.reb_integrator_trace_switch_encounter_line, TRACEKF) - elif func == "cubic": - self._S = cast(clibrebound.reb_integrator_trace_switch_encounter_predict, TRACEKF) else: self._Sfp = TRACEKF(func) self._S = self._Sfp diff --git a/src/rebound.h b/src/rebound.h index 2d802b8e8..d885e3f35 100644 --- a/src/rebound.h +++ b/src/rebound.h @@ -853,8 +853,9 @@ DLLEXPORT int reb_integrator_trace_switch_peri_fdot(struct reb_simulation* const 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); DLLEXPORT int reb_integrator_trace_switch_default(struct reb_simulation* const r, const unsigned int i, const unsigned int j); -DLLEXPORT int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r, const unsigned int i, const unsigned int j); +//DLLEXPORT int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r, const unsigned int i, const unsigned int j); DLLEXPORT int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, const unsigned int i, const unsigned int j); +DLLEXPORT int reb_integrator_trace_switch_encounter_line_peri(struct reb_simulation* const r, const unsigned int j); // debugging only // Built in collision resolve functions From 74562ebfe80eed14c03270b51053a86bb80d1cd4 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Sat, 27 Jul 2024 10:30:50 -0400 Subject: [PATCH 09/16] remove encounter predict --- src/integrator_trace.c | 248 +++++++++++++---------------------------- src/rebound.h | 1 - 2 files changed, 76 insertions(+), 173 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index 6bd3a0fb1..57916ea96 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -75,145 +75,23 @@ int reb_integrator_trace_switch_default(struct reb_simulation* const r, const un return d2*d2*d2 < dcritmax6; } -int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r, const unsigned int i, const unsigned int j){ - const struct reb_integrator_trace* const ri_trace = &(r->ri_trace); - const double dt = r->dt; - const double h2 = dt/2.; +int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, const unsigned int i, const unsigned int j){ + struct reb_integrator_trace* const ri_trace = &(r->ri_trace); + const double h2 = r->dt/2.; - double GM = r->G*r->particles[0].m; // Think this is still the correct mass, for these encounters solar mass should be more relevant - - // current - const double dxi = r->particles[i].x; // in dh + const double dxi = r->particles[i].x; const double dyi = r->particles[i].y; const double dzi = r->particles[i].z; - const double dxj = r->particles[j].x; // in dh + + const double dxj = r->particles[j].x; const double dyj = r->particles[j].y; const double dzj = r->particles[j].z; - const double dxp = dxi - dxj; - const double dyp = dyi - dyj; - const double dzp = dzi - dzj; - const double rp = (dxp*dxp + dyp*dyp + dzp*dzp); - - const double di2 = dxi*dxi + dyi*dyi + dzi*dzi; - const double di = sqrt(di2); - const double dj2 = dxj*dxj + dyj*dyj + dzj*dzj; - const double dj = sqrt(dj2); - - // second derivative - double prefacti2 = -GM/(di2*di); - double ddxi = prefacti2*dxi; - double ddyi = prefacti2*dyi; - double ddzi = prefacti2*dzi; - - double prefactj2 = -GM/(dj2*dj); - double ddxj = prefactj2*dxj; - double ddyj = prefactj2*dyj; - double ddzj = prefactj2*dzj; - - // TLUEPP - - const double pi_pre_x = r->particles[i].x - h2 * r->particles[i].vx + 0.5 * h2 * h2 * ddxi; - const double pi_pre_y = r->particles[i].y - h2 * r->particles[i].vy + 0.5 * h2 * h2 * ddyi; - const double pi_pre_z = r->particles[i].z - h2 * r->particles[i].vz + 0.5 * h2 * h2 * ddzi; - const double pi_pre_vx = r->particles[i].vx - h2 * ddxi; - const double pi_pre_vy = r->particles[i].vy - h2 * ddyi; - const double pi_pre_vz = r->particles[i].vz - h2 * ddzi; - - const double pi_post_x = r->particles[i].x + h2 * r->particles[i].vx + 0.5 * h2 * h2 * ddxi; - const double pi_post_y = r->particles[i].y + h2 * r->particles[i].vy + 0.5 * h2 * h2 * ddyi; - const double pi_post_z = r->particles[i].z + h2 * r->particles[i].vz + 0.5 * h2 * h2 * ddzi; - const double pi_post_vx = r->particles[i].vx + h2 * ddxi; - const double pi_post_vy = r->particles[i].vy + h2 * ddyi; - const double pi_post_vz = r->particles[i].vz + h2 * ddzi; - - const double pj_pre_x = r->particles[j].x - h2 * r->particles[j].vx + 0.5 * h2 * h2 * ddxj; - const double pj_pre_y = r->particles[j].y - h2 * r->particles[j].vy + 0.5 * h2 * h2 * ddyj; - const double pj_pre_z = r->particles[j].z - h2 * r->particles[j].vz + 0.5 * h2 * h2 * ddzj; - const double pj_pre_vx = r->particles[j].vx - h2 * ddxj; - const double pj_pre_vy = r->particles[j].vy - h2 * ddyj; - const double pj_pre_vz = r->particles[j].vz - h2 * ddzj; - - const double pj_post_x = r->particles[j].x + h2 * r->particles[j].vx + 0.5 * h2 * h2 * ddxj; - const double pj_post_y = r->particles[j].y + h2 * r->particles[j].vy + 0.5 * h2 * h2 * ddyj; - const double pj_post_z = r->particles[j].z + h2 * r->particles[j].vz + 0.5 * h2 * h2 * ddzj; - const double pj_post_vx = r->particles[j].vx + h2 * ddxj; - const double pj_post_vy = r->particles[j].vy + h2 * ddyj; - const double pj_post_vz = r->particles[j].vz + h2 * ddzj; - // pre - const double dx0 = pi_pre_x - pj_pre_x; - const double dy0 = pi_pre_y - pj_pre_y; - const double dz0 = pi_pre_z - pj_pre_z; - const double r0 = (dx0*dx0 + dy0*dy0 + dz0*dz0); - - // post - const double dxn = pi_post_x - pj_post_x; - const double dyn = pi_post_y - pj_post_y; - const double dzn = pi_post_z - pj_post_z; - const double rn = (dxn*dxn + dyn*dyn + dzn*dzn); - - double rmin = MIN(MIN(rn, rp), r0); - - // Quartic Polynomial - const double a = 2.*(8.*rp + dt*(rn - r0) - 4.*(rn - r0)) / (dt*dt*dt*dt); - const double b = (dt*(rn + r0) - 2.*(rn - r0)) / (dt*dt*dt); - const double c = (8.*(rn + r0) - 16.*rp - dt*(rn-r0)) / (2*dt*dt); - const double d = (6.*(rn - r0) - dt*(rn + r0)) / (4. * dt); - const double e = rp; - - // Derivative -- cubic - const double ca = 4.*a; - const double cb = 3.*b; - const double cc = 2.*c; - const double cd = d; + const double dx = dxi - dxj; + const double dy = dyi - dyj; + const double dz = dzi - dzj; + const double rp = dx*dx + dy*dy + dz*dz; - // Depressed cubic - const double p = (3.*ca*cc - cb*cb)/(3*ca*ca); - const double q = (2.*cb*cb*cb - 9.*ca*cb*cc + 27.*ca*ca*cd)/(27.*ca*ca*ca); - double roots[3] = {}; - unsigned int nroots = 0; // real roots - if (fabs(p) < 1e-10){ // p = 0 -> t^3 = -q - roots[0] = pow(-1.*q, 1./3.); - nroots = 1; - } - else if (fabs(q) < 1e-10 && p < 0){ // q = 0 -> t^3 + pt = 0 - roots[0] = 0.; - roots[1] = sqrt(-1.*p); - roots[2] = -1.*sqrt(-1.*p); - nroots = 3; - } - else{ - const double disc = q*q/4. + p*p*p/27.; - if (fabs(disc) < 1e-10){ - roots[0] = -1.5*p*q; - roots[1] = 3.*q/p; - nroots = 2; - } - else if (disc > 0.){ - const double u = pow(-q/2. - sqrt(d), 1./3.); - roots[0] = u - p/(3.*u); - nroots = 1; - } - else{ - const double u = 2.*sqrt(-p/3.); - const double t = acos(3.*q/p/u)/3.; - const double k = 2.*M_PI/3.; - roots[0] = u*cos(t); - roots[1] = u*cos(t-k); - roots[2] = u*cos(t-2*k); - nroots = 3; - } - } - - // check if roots are in the range - for (unsigned int i = 0; i < nroots; i++){ - double root = roots[i] - b/(3.*a); - if (root < h2 && root > -1.*h2){ - double rval = a*root*root*root*root + b*root*root*root + c*root*root + d*root + e; - rmin = MIN(rmin, rval); - } - } - double dcriti6 = 0.0; double dcritj6 = 0.0; @@ -231,30 +109,10 @@ int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r dcritj6 = dj2*dj2*dj2*mr*mr; } - double r_crit_hill2 = ri_trace->r_crit_hill*ri_trace->r_crit_hill; double dcritmax6 = r_crit_hill2 * r_crit_hill2 * r_crit_hill2 * MAX(dcriti6,dcritj6); - return rmin*rmin*rmin < dcritmax6; - -} - -int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, const unsigned int i, const unsigned int j){ - struct reb_integrator_trace* const ri_trace = &(r->ri_trace); - const double h2 = r->dt/2.; - - const double dxi = r->particles[i].x; - const double dyi = r->particles[i].y; - const double dzi = r->particles[i].z; - - const double dxj = r->particles[j].x; - const double dyj = r->particles[j].y; - const double dzj = r->particles[j].z; - - const double dx = dxi - dxj; - const double dy = dyi - dyj; - const double dz = dzi - dzj; - const double rp = dx*dx + dy*dy + dz*dz; + if (rp*rp*rp < dcritmax6) return 1; const double dvx = r->particles[i].vx - r->particles[j].vx; const double dvy = r->particles[i].vy - r->particles[j].vy; @@ -275,38 +133,84 @@ int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, c double rmin1 = MIN(rpre,rpost); double rmin = MIN(rmin1,rp); + if (rmin*rmin*rmin < dcritmax6) return 1; + const double a = dvx*dvx + dvy*dvy + dvz*dvz; const double b = 2.*(dx*dvx+dy*dvy+dz*dvz); const double c = dx*dx + dy*dy + dz*dz; - const double tmin = -b/(2*a); - if (tmin > -1.*h2 && tmin < h2){ - double rmin2 = a*tmin*tmin + b*tmin + c; + const double tmin_plus = -b/(2*a); + const double tmin_minus= b/(2*a); + + if (tmin_plus < h2 && tmin_plus > 0.){ + double rmin2 = a*tmin_plus*tmin_plus + b*tmin_plus + c; rmin = MIN(rmin2, rmin); } + if (tmin_minus > -1.*h2 && tmin_minus < 0.){ + double rmin2 = a*tmin_minus*tmin_minus - b*tmin_minus + c; + rmin = MIN(rmin2, rmin); + } - double dcriti6 = 0.0; - double dcritj6 = 0.0; - - const double m0 = r->particles[0].m; + return rmin*rmin*rmin < dcritmax6; +} - if (r->particles[i].m != 0){ - const double di2 = dxi*dxi + dyi*dyi + dzi*dzi; - const double mr = r->particles[i].m/(3.*m0); - dcriti6 = di2*di2*di2*mr*mr; +// Debugging only +int reb_integrator_trace_switch_encounter_line_peri(struct reb_simulation* const r, const unsigned int j){ + struct reb_integrator_trace* const ri_trace = &(r->ri_trace); + const double peri = ri_trace->peri_crit_distance; + if (peri == 0.0){ + reb_simulation_warning(r,"Pericenter distance parameter not set. Set with r->ri_trace.peri_crit_distance="); } + + const double h2 = r->dt/2.; - if (r->particles[j].m != 0){ - const double dj2 = dxj*dxj + dyj*dyj + dzj*dzj; - const double mr = r->particles[j].m/(3.*m0); - dcritj6 = dj2*dj2*dj2*mr*mr; - } + const double dx = r->particles[0].x - r->particles[j].x; + const double dy = r->particles[0].y - r->particles[j].y; + const double dz = r->particles[0].z - r->particles[j].z; + const double d2 = dx*dx + dy*dy + dz*dz; - double r_crit_hill2 = ri_trace->r_crit_hill*ri_trace->r_crit_hill; - double dcritmax6 = r_crit_hill2 * r_crit_hill2 * r_crit_hill2 * MAX(dcriti6,dcritj6); + if (d2 < peri*peri) return 1; + + const double vx = r->particles[j].vx; + const double vy = r->particles[j].vy; + const double vz = r->particles[j].vz; - return rmin*rmin*rmin < dcritmax6; + // pre + const double dx_pre = dx - h2 * vx; + const double dy_pre = dy - h2 * vy; + const double dz_pre = dz - h2 * vz; + const double dpre2 = dx_pre*dx_pre + dy_pre*dy_pre + dz_pre*dz_pre; + + // post + const double dx_post = dx + h2 * vx; + const double dy_post = dy + h2 * vy; + const double dz_post = dz + h2 * vz; + const double dpost2 = dx_post*dx_post + dy_post*dy_post + dz_post*dz_post; + + double rmin1 = MIN(dpre2,dpost2); + double rmin = MIN(rmin1,d2); + + if (rmin < peri*peri) return 1; + + const double a = vx*vx + vy*vy + vz*vz; + const double b = 2.*(dx*vx+dy*vy+dz*vz); + const double c = dx*dx + dy*dy + dz*dz; + + const double tmin_plus = -b/(2*a); + const double tmin_minus= b/(2*a); + + if (tmin_plus < h2 && tmin_plus > 0.){ + double rmin2 = a*tmin_plus*tmin_plus + b*tmin_plus + c; + rmin = MIN(rmin2, rmin); + } + + if (tmin_minus > -1.*h2 && tmin_minus < 0.){ + double rmin2 = a*tmin_minus*tmin_minus - b*tmin_minus + c; + rmin = MIN(rmin2, rmin); + } + + return rmin < peri*peri; } int reb_integrator_trace_switch_peri_distance(struct reb_simulation* const r, const unsigned int j){ diff --git a/src/rebound.h b/src/rebound.h index d885e3f35..f055369a9 100644 --- a/src/rebound.h +++ b/src/rebound.h @@ -853,7 +853,6 @@ DLLEXPORT int reb_integrator_trace_switch_peri_fdot(struct reb_simulation* const 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); DLLEXPORT int reb_integrator_trace_switch_default(struct reb_simulation* const r, const unsigned int i, const unsigned int j); -//DLLEXPORT int reb_integrator_trace_switch_encounter_predict(struct reb_simulation* const r, const unsigned int i, const unsigned int j); DLLEXPORT int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, const unsigned int i, const unsigned int j); DLLEXPORT int reb_integrator_trace_switch_encounter_line_peri(struct reb_simulation* const r, const unsigned int j); // debugging only From 0a6866756c7f271000237110b7face4b057b21f1 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Sun, 28 Jul 2024 19:17:08 -0400 Subject: [PATCH 10/16] try hd24 function --- src/integrator_trace.c | 59 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 57 insertions(+), 2 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index 57916ea96..f7979398b 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -132,9 +132,64 @@ int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, c double rmin1 = MIN(rpre,rpost); double rmin = MIN(rmin1,rp); - + if (rmin*rmin*rmin < dcritmax6) return 1; + const double dpx = r->particles[i].m*r->particles[i].vx - r->particles[j].m*r->particles[j].vx; + const double dpy = r->particles[i].m*r->particles[i].vy - r->particles[j].m*r->particles[j].vy; + const double dpz = r->particles[i].m*r->particles[i].vz - r->particles[j].m*r->particles[j].vz; + const double p2 = dpx*dpx + dpy*dpy + dpz*dpz; + + const double qp = dx*dpx + dy*dpy + dz*dpz; + double tmin; + + if (qp < 0){ + // positive solution + tmin = -1.*qp / p2; + if (tmin < h2 && tmin > 0.){ + const double rmin2 = rp + tmin; + rmin = MIN(rmin2, rmin); + } + else if (tmin > h2){ + const double rmin2 = rp + 2.*r->dt*qp + p2 + r->dt*r->dt; + rmin = MIN(rmin2, rmin); + } + } + else{ + // negative solution + tmin = qp / (dpx*dpx + dpy*dpy + dpz*dpz); + if (tmin > -1.*h2 && tmin < 0.){ + const double rmin2 = rp + tmin; + rmin = MIN(rmin2, rmin); + } + else if (tmin < -1.*h2){ + const double rmin2 = rp - 2.*r->dt*qp + p2 + r->dt*r->dt; + rmin = MIN(rmin2, rmin); + } + } + + +/* + const double a = 2.*(rpost - 2.*rpre + rp)/(r->dt*r->dt); + const double a_min = 2.*(rpre - 2.*rpost + rp)/(r->dt*r->dt); + const double b = (rpre-rpost)/(r->dt); + const double b_min = (rpost-rpre)/(r->dt); + const double c = rp; + + const double tmin_plus = -b/(2*a); + const double tmin_min = -b_min/(2*a_min); + + if (tmin_plus > 0. && tmin_plus < h2){ + double rmin2 = a*tmin_plus*tmin_plus + b*tmin_plus + c; + rmin = MIN(rmin2, rmin); + } + + if (tmin_min > -1.*h2 && tmin_min < 0.){ + double rmin2 = a*tmin_min*tmin_min + b*tmin_min + c; + rmin = MIN(rmin2, rmin); + } +*/ + /* const double a = dvx*dvx + dvy*dvy + dvz*dvz; const double b = 2.*(dx*dvx+dy*dvy+dz*dvz); const double c = dx*dx + dy*dy + dz*dz; @@ -151,7 +206,7 @@ int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, c double rmin2 = a*tmin_minus*tmin_minus - b*tmin_minus + c; rmin = MIN(rmin2, rmin); } - +*/ return rmin*rmin*rmin < dcritmax6; } From 41744a5636e966f366d4dbb2bb3d218f42022cb5 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Wed, 31 Jul 2024 16:08:18 -0400 Subject: [PATCH 11/16] correcting typos --- src/integrator_trace.c | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index f7979398b..a0f6042cc 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -145,25 +145,25 @@ int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, c if (qp < 0){ // positive solution - tmin = -1.*qp / p2; + tmin = -1.*qp/p2; if (tmin < h2 && tmin > 0.){ - const double rmin2 = rp + tmin; + const double rmin2 = rp - qp*qp/p2; rmin = MIN(rmin2, rmin); } else if (tmin > h2){ - const double rmin2 = rp + 2.*r->dt*qp + p2 + r->dt*r->dt; + const double rmin2 = rp + 2.*r->dt*qp + p2*r->dt*r->dt; rmin = MIN(rmin2, rmin); } } else{ // negative solution - tmin = qp / (dpx*dpx + dpy*dpy + dpz*dpz); + tmin = qp/p2; if (tmin > -1.*h2 && tmin < 0.){ - const double rmin2 = rp + tmin; + const double rmin2 = rp - qp*qp/p2; rmin = MIN(rmin2, rmin); } else if (tmin < -1.*h2){ - const double rmin2 = rp - 2.*r->dt*qp + p2 + r->dt*r->dt; + const double rmin2 = rp - 2.*r->dt*qp + p2*r->dt*r->dt; rmin = MIN(rmin2, rmin); } } From 9ed9bf1426c46cc34d8b627d7ce18fd2744ef372 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Wed, 31 Jul 2024 18:01:59 -0400 Subject: [PATCH 12/16] Make LINE the default switching conditions + minor addition to docs --- docs/index.md | 2 +- docs/integrators.md | 2 +- rebound/integrators/trace.py | 2 - rebound/tests/test_trace_basic.py | 3 +- src/integrator_trace.c | 93 ------------------------------- src/rebound.h | 1 - 6 files changed, 3 insertions(+), 100 deletions(-) diff --git a/docs/index.md b/docs/index.md index 154f5656e..f0e6490c2 100644 --- a/docs/index.md +++ b/docs/index.md @@ -10,8 +10,8 @@ REBOUND is an N-body integrator, i.e. a software package that can integrate the * No dependencies on external libraries. * Runs natively on Linux, MacOS, and Windows. * Symplectic integrators ([WHFast](integrators/#whfast), [SEI](integrators/#sei), [LEAPFROG](integrators/#leapfrog), [EOS](integrators/#embedded-operator-splitting-method-eos)) -* Hybrid symplectic integrators for planetary dynamics with close encounters ([MERCURIUS](integrators/#mercurius)) * Hybrid reversible integrators for planetary dynamics with arbitrary close encounters ([TRACE](integrators/#trace)) +* Hybrid symplectic integrators for planetary dynamics with close encounters ([MERCURIUS](integrators/#mercurius)) * High order symplectic integrators for integrating planetary systems ([SABA](integrators/#saba), WH Kernel methods) * High accuracy non-symplectic integrator with adaptive time-stepping ([IAS15](integrators/#ias15)) * Can integrate arbitrary user-defined ODEs that are coupled to N-body dynamics for tides, spin, etc diff --git a/docs/integrators.md b/docs/integrators.md index bfabb25ac..8d97a5b74 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. -A paper describing the TRACE implementation is in preparation. +The TRACE implementation is described in [Lu, Hernandez & Rein](https://ui.adsabs.harvard.edu/abs/2024arXiv240503800L/abstract). The following code enables TRACE and sets the critical radius to 4 Hill radii diff --git a/rebound/integrators/trace.py b/rebound/integrators/trace.py index a78acb840..e42d746db 100644 --- a/rebound/integrators/trace.py +++ b/rebound/integrators/trace.py @@ -74,8 +74,6 @@ def S(self): def S(self, func): if func == "default": self._S = cast(clibrebound.reb_integrator_trace_switch_default,TRACEKF) - elif func == "line": - self._S = cast(clibrebound.reb_integrator_trace_switch_encounter_line, TRACEKF) else: self._Sfp = TRACEKF(func) self._S = self._Sfp diff --git a/rebound/tests/test_trace_basic.py b/rebound/tests/test_trace_basic.py index ec8991a5a..8691e9361 100644 --- a/rebound/tests/test_trace_basic.py +++ b/rebound/tests/test_trace_basic.py @@ -40,11 +40,10 @@ def test_trace_encounter_prediction(self): sim = rebound.Simulation() sim.add(m=1) sim.add(m=9.55e-4,x=5.2) - sim.add(x=5.3,y=0.36,vy=-7.2) # Normal switching condition misses this + sim.add(x=5.3,y=0.36,vy=-7.2) # Non-encounter prediction misses this sim.integrator = "TRACE" sim.dt = 0.01 sim.ri_trace.r_crit_hill = 1 - sim.S = "line" sim.step() self.assertEqual(sim.ri_trace._encounter_N,3) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index a0f6042cc..bde951a9a 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -40,42 +40,6 @@ #define MAX(a, b) ((a) > (b) ? (a) : (b)) ///< Returns the maximum of a and b int reb_integrator_trace_switch_default(struct reb_simulation* const r, const unsigned int i, const unsigned int j){ - struct reb_integrator_trace* const ri_trace = &(r->ri_trace); - double dcriti6 = 0.0; - double dcritj6 = 0.0; - - const double m0 = r->particles[0].m; - - if (r->particles[i].m != 0){ - const double dxi = r->particles[i].x; // in dh - const double dyi = r->particles[i].y; - const double dzi = r->particles[i].z; - const double di2 = dxi*dxi + dyi*dyi + dzi*dzi; - const double mr = r->particles[i].m/(3.*m0); - dcriti6 = di2*di2*di2*mr*mr; - } - - if (r->particles[j].m != 0){ - const double dxj = r->particles[j].x; // in dh - const double dyj = r->particles[j].y; - const double dzj = r->particles[j].z; - const double dj2 = dxj*dxj + dyj*dyj + dzj*dzj; - const double mr = r->particles[j].m/(3.*m0); - dcritj6 = dj2*dj2*dj2*mr*mr; - } - - const double dx = r->particles[i].x - r->particles[j].x; - const double dy = r->particles[i].y - r->particles[j].y; - const double dz = r->particles[i].z - r->particles[j].z; - const double d2 = dx*dx + dy*dy + dz*dz; - - double r_crit_hill2 = ri_trace->r_crit_hill*ri_trace->r_crit_hill; - double dcritmax6 = r_crit_hill2 * r_crit_hill2 * r_crit_hill2 * MAX(dcriti6,dcritj6); - - return d2*d2*d2 < dcritmax6; -} - -int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, const unsigned int i, const unsigned int j){ struct reb_integrator_trace* const ri_trace = &(r->ri_trace); const double h2 = r->dt/2.; @@ -210,63 +174,6 @@ int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, c return rmin*rmin*rmin < dcritmax6; } -// Debugging only -int reb_integrator_trace_switch_encounter_line_peri(struct reb_simulation* const r, const unsigned int j){ - struct reb_integrator_trace* const ri_trace = &(r->ri_trace); - const double peri = ri_trace->peri_crit_distance; - if (peri == 0.0){ - reb_simulation_warning(r,"Pericenter distance parameter not set. Set with r->ri_trace.peri_crit_distance="); - } - - const double h2 = r->dt/2.; - - const double dx = r->particles[0].x - r->particles[j].x; - const double dy = r->particles[0].y - r->particles[j].y; - const double dz = r->particles[0].z - r->particles[j].z; - const double d2 = dx*dx + dy*dy + dz*dz; - - if (d2 < peri*peri) return 1; - - const double vx = r->particles[j].vx; - const double vy = r->particles[j].vy; - const double vz = r->particles[j].vz; - - // pre - const double dx_pre = dx - h2 * vx; - const double dy_pre = dy - h2 * vy; - const double dz_pre = dz - h2 * vz; - const double dpre2 = dx_pre*dx_pre + dy_pre*dy_pre + dz_pre*dz_pre; - - // post - const double dx_post = dx + h2 * vx; - const double dy_post = dy + h2 * vy; - const double dz_post = dz + h2 * vz; - const double dpost2 = dx_post*dx_post + dy_post*dy_post + dz_post*dz_post; - - double rmin1 = MIN(dpre2,dpost2); - double rmin = MIN(rmin1,d2); - - if (rmin < peri*peri) return 1; - - const double a = vx*vx + vy*vy + vz*vz; - const double b = 2.*(dx*vx+dy*vy+dz*vz); - const double c = dx*dx + dy*dy + dz*dz; - - const double tmin_plus = -b/(2*a); - const double tmin_minus= b/(2*a); - - if (tmin_plus < h2 && tmin_plus > 0.){ - double rmin2 = a*tmin_plus*tmin_plus + b*tmin_plus + c; - rmin = MIN(rmin2, rmin); - } - - if (tmin_minus > -1.*h2 && tmin_minus < 0.){ - double rmin2 = a*tmin_minus*tmin_minus - b*tmin_minus + c; - rmin = MIN(rmin2, rmin); - } - - return rmin < peri*peri; -} int reb_integrator_trace_switch_peri_distance(struct reb_simulation* const r, const unsigned int j){ const struct reb_integrator_trace* const ri_trace = &(r->ri_trace); diff --git a/src/rebound.h b/src/rebound.h index f055369a9..88f41c7e7 100644 --- a/src/rebound.h +++ b/src/rebound.h @@ -853,7 +853,6 @@ DLLEXPORT int reb_integrator_trace_switch_peri_fdot(struct reb_simulation* const 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); DLLEXPORT int reb_integrator_trace_switch_default(struct reb_simulation* const r, const unsigned int i, const unsigned int j); -DLLEXPORT int reb_integrator_trace_switch_encounter_line(struct reb_simulation* const r, const unsigned int i, const unsigned int j); DLLEXPORT int reb_integrator_trace_switch_encounter_line_peri(struct reb_simulation* const r, const unsigned int j); // debugging only From 2c9841336ad86aa109bc591965ea6c6bad57dd35 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Sun, 11 Aug 2024 13:08:56 +0200 Subject: [PATCH 13/16] typos in switching condition --- src/integrator_trace.c | 104 ++++++++--------------------------------- 1 file changed, 20 insertions(+), 84 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index bde951a9a..e3abcaf45 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -81,97 +81,33 @@ int reb_integrator_trace_switch_default(struct reb_simulation* const r, const un const double dvx = r->particles[i].vx - r->particles[j].vx; const double dvy = r->particles[i].vy - r->particles[j].vy; const double dvz = r->particles[i].vz - r->particles[j].vz; - - // pre - const double dx_pre = dx - h2 * dvx; - const double dy_pre = dy - h2 * dvy; - const double dz_pre = dz - h2 * dvz; - const double rpre = dx_pre*dx_pre + dy_pre*dy_pre + dz_pre*dz_pre; - - // post - const double dx_post = dx + h2 * dvx; - const double dy_post = dy + h2 * dvy; - const double dz_post = dz + h2 * dvz; - const double rpost = dx_post*dx_post + dy_post*dy_post + dz_post*dz_post; - - double rmin1 = MIN(rpre,rpost); - double rmin = MIN(rmin1,rp); + const double v2 = dvx*dvx + dvy*dvy + dvz*dvz; - if (rmin*rmin*rmin < dcritmax6) return 1; - - const double dpx = r->particles[i].m*r->particles[i].vx - r->particles[j].m*r->particles[j].vx; - const double dpy = r->particles[i].m*r->particles[i].vy - r->particles[j].m*r->particles[j].vy; - const double dpz = r->particles[i].m*r->particles[i].vz - r->particles[j].m*r->particles[j].vz; - const double p2 = dpx*dpx + dpy*dpy + dpz*dpz; - - const double qp = dx*dpx + dy*dpy + dz*dpz; - double tmin; - - if (qp < 0){ - // positive solution - tmin = -1.*qp/p2; - if (tmin < h2 && tmin > 0.){ - const double rmin2 = rp - qp*qp/p2; - rmin = MIN(rmin2, rmin); - } - else if (tmin > h2){ - const double rmin2 = rp + 2.*r->dt*qp + p2*r->dt*r->dt; - rmin = MIN(rmin2, rmin); - } + const double qv = dx*dvx + dy*dvy + dz*dvz; + int d; + + if (qv == 0.0){ // Small + // minimum is at present, which is already checked for + return 0; + } + else if (qv < 0){ + d = 1; } else{ - // negative solution - tmin = qp/p2; - if (tmin > -1.*h2 && tmin < 0.){ - const double rmin2 = rp - qp*qp/p2; - rmin = MIN(rmin2, rmin); - } - else if (tmin < -1.*h2){ - const double rmin2 = rp - 2.*r->dt*qp + p2*r->dt*r->dt; - rmin = MIN(rmin2, rmin); - } + d = -1; } - -/* - const double a = 2.*(rpost - 2.*rpre + rp)/(r->dt*r->dt); - const double a_min = 2.*(rpre - 2.*rpost + rp)/(r->dt*r->dt); - const double b = (rpre-rpost)/(r->dt); - const double b_min = (rpost-rpre)/(r->dt); - const double c = rp; - - const double tmin_plus = -b/(2*a); - const double tmin_min = -b_min/(2*a_min); - - if (tmin_plus > 0. && tmin_plus < h2){ - double rmin2 = a*tmin_plus*tmin_plus + b*tmin_plus + c; - rmin = MIN(rmin2, rmin); + double dmin2; + double tmin = -d*qv/(2.*v2); + if (fabs(tmin) < h2){ + // minimum is in the window + dmin2 = rp - qv/(4.*v2); } - - if (tmin_min > -1.*h2 && tmin_min < 0.){ - double rmin2 = a*tmin_min*tmin_min + b*tmin_min + c; - rmin = MIN(rmin2, rmin); - } -*/ - /* - const double a = dvx*dvx + dvy*dvy + dvz*dvz; - const double b = 2.*(dx*dvx+dy*dvy+dz*dvz); - const double c = dx*dx + dy*dy + dz*dz; - - const double tmin_plus = -b/(2*a); - const double tmin_minus= b/(2*a); - - if (tmin_plus < h2 && tmin_plus > 0.){ - double rmin2 = a*tmin_plus*tmin_plus + b*tmin_plus + c; - rmin = MIN(rmin2, rmin); - } - - if (tmin_minus > -1.*h2 && tmin_minus < 0.){ - double rmin2 = a*tmin_minus*tmin_minus - b*tmin_minus + c; - rmin = MIN(rmin2, rmin); + else{ + dmin2 = rp + d*qv*h2 + v2*h2*h2/4.; } -*/ - return rmin*rmin*rmin < dcritmax6; + + return dmin2*dmin2*dmin2 < dcritmax6; } From 2c0514410ee07d5d22c6751212ad5aa7aa2024c0 Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Mon, 12 Aug 2024 09:43:50 +0200 Subject: [PATCH 14/16] removed absolute values in switching condition --- src/integrator_trace.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index e3abcaf45..c903294fc 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -99,7 +99,7 @@ int reb_integrator_trace_switch_default(struct reb_simulation* const r, const un double dmin2; double tmin = -d*qv/(2.*v2); - if (fabs(tmin) < h2){ + if (tmin < h2){ // minimum is in the window dmin2 = rp - qv/(4.*v2); } From 15166627fb4d49eac720ba2a8284bd4080c387ed Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Wed, 14 Aug 2024 10:39:20 +0200 Subject: [PATCH 15/16] factor of 2 --- src/integrator_trace.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/integrator_trace.c b/src/integrator_trace.c index c903294fc..a93f856cf 100644 --- a/src/integrator_trace.c +++ b/src/integrator_trace.c @@ -98,13 +98,13 @@ int reb_integrator_trace_switch_default(struct reb_simulation* const r, const un } double dmin2; - double tmin = -d*qv/(2.*v2); + double tmin = -d*qv/v2; if (tmin < h2){ // minimum is in the window - dmin2 = rp - qv/(4.*v2); + dmin2 = rp - qv*qv/v2; } else{ - dmin2 = rp + d*qv*h2 + v2*h2*h2/4.; + dmin2 = rp + 2*d*qv*h2 + v2*h2*h2; } return dmin2*dmin2*dmin2 < dcritmax6; From 7f641fb794f6fad481371f27fdce68a39e03447d Mon Sep 17 00:00:00 2001 From: Tiger Lu Date: Thu, 15 Aug 2024 15:35:47 +0200 Subject: [PATCH 16/16] added back default pericenter condition --- src/rebound.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rebound.h b/src/rebound.h index 88f41c7e7..e1f239697 100644 --- a/src/rebound.h +++ b/src/rebound.h @@ -849,11 +849,11 @@ DLLEXPORT double reb_integrator_mercurius_L_C5(const struct reb_simulation* cons // Built in trace switching functions +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); DLLEXPORT int reb_integrator_trace_switch_default(struct reb_simulation* const r, const unsigned int i, const unsigned int j); -DLLEXPORT int reb_integrator_trace_switch_encounter_line_peri(struct reb_simulation* const r, const unsigned int j); // debugging only // Built in collision resolve functions