Skip to content

Commit 4f912df

Browse files
hannoreinhannorein
committed
ias
1 parent 7ce2098 commit 4f912df

1 file changed

Lines changed: 52 additions & 16 deletions

File tree

‎src/integrator_ias15.c‎

Lines changed: 52 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -170,7 +170,7 @@ void reb_integrator_ias15_alloc(struct reb_simulation* r){
170170
}else if (r->integrator==REB_INTEGRATOR_TRACE && r->ri_trace.mode == REB_TRACE_MODE_KEPLER){
171171
N3 = 3*r->ri_trace.encounter_N;// trace close encounter
172172
}else{
173-
N3 = 3*r->N;
173+
N3 = 3*(r->N+r->N_varX);
174174
}
175175
if (N3 > r->ri_ias15.N_allocated) {
176176
realloc_dp7(&(r->ri_ias15.g),N3);
@@ -195,12 +195,12 @@ void reb_integrator_ias15_alloc(struct reb_simulation* r){
195195
}
196196
r->ri_ias15.N_allocated = N3;
197197
}
198-
if (N3/3 > r->ri_ias15.N_allocated_map){
199-
r->ri_ias15.map = realloc(r->ri_ias15.map,sizeof(size_t)*(N3/3));
200-
for (size_t i=0;i<N3/3;i++){
198+
if (r->N > r->ri_ias15.N_allocated_map){
199+
r->ri_ias15.map = realloc(r->ri_ias15.map,sizeof(size_t)*r->N);
200+
for (size_t i=0;i<r->N;i++){
201201
r->ri_ias15.map[i] = i;
202202
}
203-
r->ri_ias15.N_allocated_map = N3/3;
203+
r->ri_ias15.N_allocated_map = r->N;
204204
}
205205

206206
}
@@ -212,7 +212,9 @@ static int reb_integrator_ias15_step_try(struct reb_simulation* r) {
212212
reb_integrator_ias15_alloc(r);
213213

214214
struct reb_particle* const particles = r->particles;
215-
size_t N;
215+
struct reb_particle* const particles_varX = r->particles_varX;
216+
size_t N; // Number of real particles
217+
size_t N_varX = 0; // Number of variational particles
216218
size_t* map; // this map allow for integrating only a selection of particles
217219
if (r->integrator==REB_INTEGRATOR_MERCURIUS){// mercurius close encounter
218220
N = r->ri_mercurius.encounter_N;
@@ -230,9 +232,10 @@ static int reb_integrator_ias15_step_try(struct reb_simulation* r) {
230232
}
231233
}else{
232234
N = r->N;
235+
N_varX = r->N_varX;
233236
map = r->ri_ias15.map; // identity map
234237
}
235-
const size_t N3 = 3*N;
238+
const size_t N3 = 3*(N+N_varX);
236239

237240
reb_simulation_update_acceleration(r);
238241

@@ -262,8 +265,20 @@ static int reb_integrator_ias15_step_try(struct reb_simulation* r) {
262265
a0[3*k+1] = particles[mk].ay;
263266
a0[3*k+2] = particles[mk].az;
264267
}
268+
for(size_t mk=0;mk<N_varX;mk++) {
269+
size_t k = mk + N;
270+
x0[3*k] = particles_varX[mk].x;
271+
x0[3*k+1] = particles_varX[mk].y;
272+
x0[3*k+2] = particles_varX[mk].z;
273+
v0[3*k] = particles_varX[mk].vx;
274+
v0[3*k+1] = particles_varX[mk].vy;
275+
v0[3*k+2] = particles_varX[mk].vz;
276+
a0[3*k] = particles_varX[mk].ax;
277+
a0[3*k+1] = particles_varX[mk].ay;
278+
a0[3*k+2] = particles_varX[mk].az;
279+
}
265280
if (r->gravity==REB_GRAVITY_COMPENSATED){
266-
for(size_t k=0;k<N;k++) {
281+
for(size_t k=0;k<N+N_varX;k++) {
267282
size_t mk = map[k];
268283
csa0[3*k] = gravity_cs[mk].x;
269284
csa0[3*k+1] = gravity_cs[mk].y;
@@ -336,7 +351,7 @@ static int reb_integrator_ias15_step_try(struct reb_simulation* r) {
336351
r->t = t_beginning + r->dt * h[n];
337352

338353
// Prepare particles arrays for force calculation
339-
for(size_t i=0;i<N;i++) { // Predict positions at interval n using b values
354+
for(size_t i=0;i<N+N_varX;i++) { // Predict positions at interval n using b values
340355
size_t mi = map[i];
341356
const size_t k0 = 3*i+0;
342357
const size_t k1 = 3*i+1;
@@ -353,7 +368,7 @@ static int reb_integrator_ias15_step_try(struct reb_simulation* r) {
353368
particles[mi].z = xk2 + x0[k2];
354369
}
355370
if (r->calculate_megno || (r->additional_forces && r->force_is_velocity_dependent)){
356-
for(size_t i=0;i<N;i++) { // Predict velocities at interval n using b values
371+
for(size_t i=0;i<N+N_varX;i++) { // Predict velocities at interval n using b values
357372
size_t mi = map[i];
358373
const size_t k0 = 3*i+0;
359374
const size_t k1 = 3*i+1;
@@ -377,7 +392,7 @@ static int reb_integrator_ias15_step_try(struct reb_simulation* r) {
377392
integrator_megno_thisdt += w[n] * (r->t-r->megno_initial_t) * reb_tools_megno_deltad_delta(r);
378393
}
379394

380-
for(size_t k=0;k<N;++k) {
395+
for(size_t k=0;k<N+N_varX;++k) {
381396
size_t mk = map[k];
382397
at[3*k] = particles[mk].ax;
383398
at[3*k+1] = particles[mk].ay;
@@ -530,15 +545,12 @@ static int reb_integrator_ias15_step_try(struct reb_simulation* r) {
530545
// r->ri_ias15.adaptive_mode==REB_IAS15_PRS23
531546
// Here, the acceleration, jerk and snap are used to estimate the new timestep.
532547
// The method is described in detail in Pham, Rein, Spiegel 2023
533-
size_t Nreal = N;
534548
if (r->ri_ias15.adaptive_mode==REB_IAS15_INDIVIDUAL || r->ri_ias15.adaptive_mode==REB_IAS15_GLOBAL){ // Old adaptive timestepping methods
535549
double integrator_error = 0.0; // Try to estimate integrator error based on last polynomial
536550
if (r->ri_ias15.adaptive_mode==REB_IAS15_GLOBAL){
537551
double maxa = 0.0;
538552
double maxj = 0.0;
539-
for(size_t i=0;i<Nreal;i++){ // Looping over all particles and all 3 components of the acceleration.
540-
// Note: Before December 2020, N-N_var, was simply N. This change should make timestep choices during
541-
// close encounters more stable if variational particles are present.
553+
for(size_t i=0;i<N;i++){ // Looping over all particles and all 3 components of the acceleration.
542554
size_t mi = map[i];
543555
const double v2 = particles[mi].vx*particles[mi].vx+particles[mi].vy*particles[mi].vy+particles[mi].vz*particles[mi].vz;
544556
const double x2 = particles[mi].x*particles[mi].x+particles[mi].y*particles[mi].y+particles[mi].z*particles[mi].z;
@@ -574,7 +586,7 @@ static int reb_integrator_ias15_step_try(struct reb_simulation* r) {
574586
};
575587
}else{ // adaptive_mode >= 2 (New adaptive timestepping method, default since January 2024)
576588
double min_timescale2 = INFINITY; // note factor of dt_done**2 not included
577-
for(size_t i=0;i<Nreal;i++){
589+
for(size_t i=0;i<N;i++){
578590
double a0i = 0; // (acceleration at beginning of timestep)^2
579591
double y2 = 0; // (accaleration at end of timestep)^2
580592
double y3 = 0; // (jerk * dt_done)^2
@@ -633,6 +645,20 @@ static int reb_integrator_ias15_step_try(struct reb_simulation* r) {
633645
particles[mk].ay = a0[3*k+1];
634646
particles[mk].az = a0[3*k+2];
635647
}
648+
for(size_t mk=0;mk<N_varX;++mk) {
649+
size_t k = mk + N;
650+
particles_varX[mk].x = x0[3*k+0]; // Set initial position
651+
particles_varX[mk].y = x0[3*k+1];
652+
particles_varX[mk].z = x0[3*k+2];
653+
654+
particles_varX[mk].vx = v0[3*k+0]; // Set initial velocity
655+
particles_varX[mk].vy = v0[3*k+1];
656+
particles_varX[mk].vz = v0[3*k+2];
657+
658+
particles_varX[mk].ax = a0[3*k+0]; // Set initial acceleration
659+
particles_varX[mk].ay = a0[3*k+1];
660+
particles_varX[mk].az = a0[3*k+2];
661+
}
636662
r->dt = dt_new;
637663
if (r->dt_last_done!=0.){ // Do not predict next e/b values if this is the first time step.
638664
double ratio = r->dt/r->dt_last_done;
@@ -691,6 +717,16 @@ static int reb_integrator_ias15_step_try(struct reb_simulation* r) {
691717
particles[mk].vy = v0[3*k+1];
692718
particles[mk].vz = v0[3*k+2];
693719
}
720+
for(size_t mk=0;mk<N_varX;++mk) {
721+
size_t k = mk + N;
722+
particles_varX[mk].x = x0[3*k+0]; // Set final position
723+
particles_varX[mk].y = x0[3*k+1];
724+
particles_varX[mk].z = x0[3*k+2];
725+
726+
particles_varX[mk].vx = v0[3*k+0]; // Set final velocity
727+
particles_varX[mk].vy = v0[3*k+1];
728+
particles_varX[mk].vz = v0[3*k+2];
729+
}
694730
copybuffers(e,er,N3);
695731
copybuffers(b,br,N3);
696732
double ratio = r->dt/dt_done;

0 commit comments

Comments
 (0)