Skip to content

Commit

Permalink
replaced interpolation with x, v values from function
Browse files Browse the repository at this point in the history
  • Loading branch information
@johnwez1
johnwez1 committed Oct 20, 2024
1 parent dc60164 commit dd36123
Show file tree
Hide file tree
Showing 2 changed files with 24 additions and 133 deletions.
108 changes: 21 additions & 87 deletions galpy/util/wez_ias15.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,20 +222,18 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
double *v= (double *) malloc ( dim * sizeof(double) );
double *a= (double *) malloc ( dim * sizeof(double) );
double *xs= (double *) malloc ( dim * sizeof(double) ); //x substep
double *vs= (double *) malloc ( dim * sizeof(double) ); //v substep

double *Bs= (double *) malloc ( (order * dim) * sizeof(double) );
double *Es= (double *) malloc ( (order * dim) * sizeof(double) );
double *BDs= (double *) malloc ( (order * dim) * sizeof(double) );
double *Gs= (double *) malloc ( (order * dim) * sizeof(double) );
double *Fs= (double *) malloc ( ((order + 1) * dim) * sizeof(double) );

int maxsteps = nt; // as an initial value, we don't know how many steps we will take
double *x_steps = (double *) malloc ( dim * (maxsteps + 1) * sizeof(double) ); //storing the result before interpolation
double *v_steps = (double *) malloc ( dim * (maxsteps + 1) * sizeof(double) ); //storing the result before interpolation
double *t_steps = (double *) malloc ( (maxsteps + 1) * sizeof(double) ); //storing the result before interpolation
double hs;

int ii, jj, kk;
printf("DO THE IAS15 INTEGRATION...\n");

for (ii=0; ii < dim; ii++) {
*x++= *(yo+ii);
*v++= *(yo+dim+ii);
Expand All @@ -250,8 +248,8 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po

double diff_G;

//save_ias15(dim, x, v, result);
//result+= 2 * dim;
save_ias15(dim, x, v, result);
result+= 2 * dim;

*err= 0;

Expand Down Expand Up @@ -280,9 +278,7 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
#endif
//WHILE THERE IS TIME REMAINING, INTEGRATE A DYNAMIC TIMESTEP FORWARD AND SUBTRACT FROM THE TIME REMAINING
double time_remaining = nt * dt;
int steps = 0;

save_dummy_ias15(dim, x, v, 0.0, x_steps, v_steps, t_steps, steps, maxsteps);
int steps = 1;

while(time_remaining > 0) {
if ( interrupted ) {
Expand Down Expand Up @@ -328,7 +324,7 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
}

//at start of each step reset xs
for (int l=0; l < dim; l++) *(xs+l)= *(x+l);
//for (int l=0; l < dim; l++) *(xs+l)= *(x+l);

double max_delta_B6 = 0.0; //also = max delta G
double max_a = 0.0;
Expand Down Expand Up @@ -380,91 +376,36 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
//fix for inf values issue
double dt_required;
double correction_factor = seventhroot(precision_parameter / (max_B6/max_a));
//double correction_factor = pow(precision_parameter / (max_B6/max_a), 1.0/7.0);

printf("MAXA %f\n", max_a);
printf("MAXB6 %f\n", max_B6);
printf("CORRECTION FACTOR %f\n", correction_factor);

printf("OLD DT %f\n", dt_temp);
if (isnormal(correction_factor)){
dt_required = dt_temp * correction_factor;
} else{
dt_required = dt_temp;
}

printf("NEW DT %f\n", dt_required);
if(fabs(dt_temp) > fabs(dt_required)){
printf("REJECTED TRY AGAIN...\n");
//rejected, try again with dt required
dt = dt_required;
} else {
//accepted, update position/velocity and do next timestep with dt required
time_remaining -= fabs(dt); //will eventually get negative as we stepped forward the minimum of dt and time_remaining
to = to_temp;

update_position(x, x, v, dim, 1, dt_temp, Fs, Bs);
update_velocity(v, v, dim, 1, dt_temp, Fs , Bs);
next_sequence_Bs(1, Bs, Es, BDs, dim);

//save_ias15(dim, x, v, result);
//result+= 2 * dim;

steps += 1;

//estimate the function over the interval for any points in the t array
while(fabs(to) < fabs(t[steps]) && fabs(t[steps]) <= fabs(to_temp) && steps < nt){
hs = (fabs(t[steps]) - fabs(to))/(fabs(to_temp) - fabs(to));

if(steps > maxsteps){ //dynamic arrays must be increased
printf("...RESIZING ARRAY...\n");
//break;
printf("MAXSTEPS: %d\n", maxsteps);
update_position(xs, x, v, dim, hs, dt_temp, Fs, Bs);
update_velocity(vs, v, dim, hs, dt_temp, Fs, Bs);

double *x_steps_new = (double *) malloc ( dim * (2 * maxsteps + 1) * sizeof(double) );
double *v_steps_new = (double *) malloc ( dim * (2 * maxsteps + 1) * sizeof(double) );
double *t_steps_new = (double *) malloc ( (2 * maxsteps + 1) * sizeof(double) );
save_ias15(dim, xs, vs, result);
result+= 2 * dim;

for(int k = 0; k < steps; k += 1){
t_steps_new[k] = t_steps[k];
for(ii = 0; ii < dim; ii += 1){
x_steps_new[ii * (2 * maxsteps + 1) + k] = x_steps[ii * (maxsteps + 1) + k];
v_steps_new[ii * (2 * maxsteps + 1) + k] = v_steps[ii * (maxsteps + 1) + k];
}
//printf("%f\n", t_steps[k]);
//printf("%f\n", t_steps[k]);
}
maxsteps *= 2;
printf("NEW MAXSTEPS: %d\n", maxsteps);
free(x_steps);
free(v_steps);
free(t_steps);
x_steps = x_steps_new;
v_steps = v_steps_new;
t_steps = t_steps_new;
steps += 1;
}

//printf("END OF STEP: %d\n", steps);
//printf("STEP SIZE: %f\n", fabs(dt));
save_dummy_ias15(dim, x, v, t_steps[steps - 1] + fabs(dt), x_steps, v_steps, t_steps, steps, maxsteps);
//finally, update the next step
update_position(x, x, v, dim, 1, dt_temp, Fs, Bs);
update_velocity(v, v, dim, 1, dt_temp, Fs , Bs);
next_sequence_Bs(1, Bs, Es, BDs, dim);

to = to_temp;
dt = dt_required;
//dt = 0.028;
}
}

printf("We have max %d steps\n", maxsteps);
printf("We have taken: %d steps\n", steps);

double *t_steps_short = (double *) malloc ( (steps + 1) * sizeof(double) );
//t_steps = realloc(t_steps, (steps + 1) * sizeof(double));

//printf("=============\n");
for(int k = 0; k <= steps; k += 1){
t_steps_short[k] = t_steps[k];
//printf("%f\n", t_steps[k]);
//printf("%f\n", t_steps[k]);
}
//printf("-------------\n");
//run the interpolation over each stored dimension
save_ias15(dim, x_steps, v_steps, t_steps_short, maxsteps, steps, nt, fabs(init_dt), result);
// Back to default handler
#ifndef _WIN32
action.sa_handler= SIG_DFL;
Expand All @@ -475,19 +416,12 @@ void wez_ias15(void (*func)(double t, double *q, double *a, int nargs, struct po
free(v);
free(a);
free(xs);

//printf("AAA\n");
free(vs);
free(Fs);
free(Gs);
free(Bs);
free(Es);
free(BDs);
//printf("BBB\n");
free(x_steps);
free(v_steps);
free(t_steps);
free(t_steps_short);
//We're done
}

void update_velocity(double *v, double *v1, int dim, double h_n, double T, double * Fs , double * Bs){
Expand Down
49 changes: 3 additions & 46 deletions galpy/util/wez_ias15.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,53 +58,10 @@ void wez_ias15(void (*func)(double, double *, double *,
double, double,
double *,int *);

static inline void save_ias15(int dim, double *x_steps, double *v_steps, double *t_steps, int maxsteps, int steps, int nt, double dt, double *result){
static inline void save_ias15(int dim, double *qo, double *po, double *result){
int ii;
int jj;

double *x_steps_short = (double *) malloc ( (steps + 1) * sizeof(double) ); //storing the result before interpolation
double *v_steps_short = (double *) malloc ( (steps + 1) * sizeof(double) ); //storing the result before interpolation

gsl_interp_accel *acc_x = gsl_interp_accel_alloc();
gsl_spline *spline_x = gsl_spline_alloc(gsl_interp_cspline, (steps + 1));

gsl_interp_accel *acc_v = gsl_interp_accel_alloc();
gsl_spline *spline_v = gsl_spline_alloc(gsl_interp_cspline, (steps + 1));

for (ii=0; ii < dim; ii++){
memcpy(x_steps_short, &x_steps[ii * (maxsteps + 1)], (steps + 1) * sizeof(double));
memcpy(v_steps_short, &v_steps[ii * (maxsteps + 1)], (steps + 1) * sizeof(double));

gsl_spline_init (spline_x, t_steps, x_steps_short, steps + 1);
gsl_spline_init (spline_v, t_steps, v_steps_short, steps + 1);

for (jj = 0; jj <= nt; jj += 1){
//if (ii == 0){
// printf("%d\n", jj);
// printf("%f\n", jj * dt);
//}
result[ii + (2 * dim * jj)] = gsl_spline_eval(spline_x, jj * dt, acc_x);
result[ii + (2 * dim * jj) + dim] = gsl_spline_eval(spline_v, jj * dt, acc_v);
}
}

gsl_spline_free (spline_x);
gsl_interp_accel_free (acc_x);
gsl_spline_free (spline_v);
gsl_interp_accel_free (acc_v);
free(x_steps_short);
free(v_steps_short);
}

static inline void save_dummy_ias15(int dim, double *x, double *v, double t, double *x_steps, double *v_steps, double *t_steps, int steps, int maxsteps){
int ii;
for (ii=0; ii < dim; ii++){
x_steps[ii * (maxsteps + 1) + steps] = *x++;
}
for (ii=0; ii < dim; ii++){
v_steps[ii * (maxsteps + 1) + steps] = *v++;
}
t_steps[steps] = t;
for (ii=0; ii < dim; ii++) *result++= *qo++;
for (ii=0; ii < dim; ii++) *result++= *po++;
}

void update_velocity(double *, double *, int, double, double, double *, double *);
Expand Down

0 comments on commit dd36123

Please sign in to comment.