yohai/phigrape
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Cleaned up corrector loop

Browse source
This commit is contained in:
Yohai Meiron 2021-12-26 17:20:20 -05:00
parent 15227a9eee
commit 8b71f4fe92
1 changed files with 38 additions and 32 deletions
Download patch file Download diff file Expand all files Collapse all files

70
phigrape.cpp
View file

@ -163,7 +163,7 @@ public:
{ {
ind_act_loc.resize(n_loc); ind_act_loc.resize(n_loc);
} }
double get_minimum_time(const std::vector<double> &t, std::vector<double> &dt) double get_minimum_time(const std::vector<double> &t, const std::vector<double> &dt)
{ {
double min_t_loc, min_t; double min_t_loc, min_t;
#ifdef ETICS #ifdef ETICS
@ -249,6 +249,29 @@ inline double aarseth_step(const double eta, const double dt, const double3 &a,
return sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); return sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
} }
inline double blockize_step(double dt, double dt_prev, double min_t, double dt_min, double dt_max)
{
double dt_new = dt_prev;
if (dt < dt_min) dt_prev = dt_min;
if ((dt < dt_prev) && (dt > dt_min)) {
int power = log(dt)/M_LN2 - 1;
dt_new = pow(2.0, power);
}
if ((dt > 2*dt_new) && (fmod(min_t, 2*dt_new) == 0) && (2*dt_new <= dt_max)) dt_new *= 2;
return dt_new;
}
inline void predictor(double min_t, const int n_act, const std::vector<int> &ind_act, const std::vector<double> &t, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double3> &a, const std::vector<double3> &adot, std::vector<double3> &x_act_new, std::vector<double3> &v_act_new)
{
for (int i=0; i<n_act; i++) {
int j_act = ind_act[i];
double dt = min_t - t[j_act];
double dt2half = 0.5*dt*dt;
double dt3over6 = (1./3.)*dt*dt2half;
x_act_new[i] = x[j_act] + v[j_act]*dt + a[j_act]*dt2half + adot[j_act]*dt3over6;
v_act_new[i] = v[j_act] + a[j_act]*dt + adot[j_act]*dt2half;
} /* i */
}
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
@ -521,6 +544,8 @@ int main(int argc, char *argv[])
/* Find the BH(s) indices in the active list */ /* Find the BH(s) indices in the active list */
int i_bh1=0, i_bh2=1; int i_bh1=0, i_bh2=1;
#ifdef ETICS #ifdef ETICS
/* Unlike with the simple active search, with GPU accelerated GRAPite
active search, the list of active indices is not sorted. */
int n_bh = config.live_smbh_count; int n_bh = config.live_smbh_count;
if (config.grapite_active_search && (n_bh>0)) { if (config.grapite_active_search && (n_bh>0)) {
int act_def_grapite_bh_count = 0; int act_def_grapite_bh_count = 0;
@ -537,14 +562,7 @@ int main(int argc, char *argv[])
#endif #endif
/* predict the active particles positions etc... on all the nodes */ /* predict the active particles positions etc... on all the nodes */
for (int i=0; i<n_act; i++) { predictor(min_t, n_act, ind_act, t, x, v, a, adot, x_act_new, v_act_new);
int j_act = ind_act[i];
double dt = min_t - t[j_act];
double dt2half = 0.5*dt*dt;
double dt3over6 = (1./3.)*dt*dt2half;
x_act_new[i] = x[j_act] + v[j_act]*dt + a[j_act]*dt2half + adot[j_act]*dt3over6;
v_act_new[i] = v[j_act] + a[j_act]*dt + adot[j_act]*dt2half;
} /* i */
/* Calculate gravity on active particles */ /* Calculate gravity on active particles */
calc_self_grav(min_t, n_act, ind_act, x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new); calc_self_grav(min_t, n_act, ind_act, x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new);
@ -562,46 +580,37 @@ int main(int argc, char *argv[])
} }
/* correct the active particles positions etc... on all the nodes */ /* correct the active particles positions etc... on all the nodes */
double min_dt = dt_max; double min_dt = dt_max; // notice that min_dt is not the same as dt_min; this one is to store the minimum timestep among currently active particles
for (int i=0; i<n_act; i++) { for (int i=0; i<n_act; i++) {
int j_act = ind_act[i]; int j_act = ind_act[i];
double dt_tmp = min_t - t[j_act]; double dt_cur = dt[j_act];
double3 a2, a3; double3 a2, a3;
calc_high_derivatives(dt_tmp, a[j_act], a_act_new[i], adot[j_act], adot_act_new[i], a2, a3); calc_high_derivatives(dt_cur, a[j_act], a_act_new[i], adot[j_act], adot_act_new[i], a2, a3);
corrector(dt_tmp, a2, a3, x_act_new[i], v_act_new[i]); corrector(dt_cur, a2, a3, x_act_new[i], v_act_new[i]);
//TODO make beautiful //TODO make beautiful
double eta_curr; double eta_curr;
if ((config.live_smbh_count > 0) && (ind_act[i] < config.live_smbh_count)) eta_curr = config.eta/config.eta_bh_corr; if ((config.live_smbh_count > 0) && (ind_act[i] < config.live_smbh_count)) eta_curr = config.eta/config.eta_bh_corr;
else eta_curr = config.eta; else eta_curr = config.eta;
double dt_new = aarseth_step(eta_curr, dt_tmp, a_act_new[i], adot_act_new[i], a2, a3); double dt_new = aarseth_step(eta_curr, dt_cur, a_act_new[i], adot_act_new[i], a2, a3);
//TODO the below should be moved to a function dt_new = blockize_step(dt_new, dt_cur, min_t, dt_min, dt_max);
if (dt_new < dt_min) dt_tmp = dt_min;
if ((dt_new < dt_tmp) && (dt_new > dt_min)) {
int power = log(dt_new)/log(2.0) - 1;
dt_tmp = pow(2.0, (double)power); // TODO why is this casting needed here?
}
if ((dt_new > 2*dt_tmp) && (fmod(min_t, 2*dt_tmp) == 0) && (2*dt_tmp <= dt_max)) {
dt_tmp *= 2;
}
if (config.dt_min_warning && (myRank == 0)) { if (config.dt_min_warning && (myRank == 0)) {
if (dt_tmp == dt_min) { if (dt_new == dt_min) {
printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d time_cur=%.16E\n", dt_tmp, ind_act[i], time_cur); printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d time_cur=%.16E\n", dt_cur, ind_act[i], time_cur);
fflush(stdout); fflush(stdout);
} }
} }
if (dt_tmp < min_dt) min_dt = dt_tmp; if (dt_new < min_dt) min_dt = dt_new;
x[j_act] = x_act_new[i]; x[j_act] = x_act_new[i];
v[j_act] = v_act_new[i]; v[j_act] = v_act_new[i];
t[j_act] = min_t; t[j_act] = min_t;
dt[j_act] = dt_tmp; dt[j_act] = dt_new;
pot[j_act] = pot_act_new[i]; pot[j_act] = pot_act_new[i];
pot_ext[j_act] = pot_act_ext[i]; pot_ext[j_act] = pot_act_ext[i];
a[j_act] = a_act_new[i]; a[j_act] = a_act_new[i];
@ -614,10 +623,8 @@ int main(int argc, char *argv[])
if (config.live_smbh_count==2) dt[1] = min_dt; if (config.live_smbh_count==2) dt[1] = min_dt;
} }
if (config.binary_smbh_influence_sphere_output && (myRank == rootRank)) { if (config.binary_smbh_influence_sphere_output && (myRank == rootRank))
//TODO clean up this mass. We don't want to have all these _act arrays; they are only needed for this lame function.
binary_smbh_influence_sphere_output(ind_act, n_act, timesteps, time_cur); binary_smbh_influence_sphere_output(ind_act, n_act, timesteps, time_cur);
}
/* load the new values for active particles to the local GRAPE's */ /* load the new values for active particles to the local GRAPE's */
for (int i=0; i<n_act; i++) { for (int i=0; i<n_act; i++) {
@ -629,7 +636,6 @@ int main(int argc, char *argv[])
int j_act = ind_act[i]; int j_act = ind_act[i];
int address = ind_act[i] - myRank*n_loc; int address = ind_act[i] - myRank*n_loc;
g6_set_j_particle(clusterid, address, ind_act[i], t[j_act], dt[j_act], m[j_act], zeros, adot[j_act]*(1./6.), a[j_act]*0.5, v[j_act], x[j_act]); g6_set_j_particle(clusterid, address, ind_act[i], t[j_act], dt[j_act], m[j_act], zeros, adot[j_act]*(1./6.), a[j_act]*0.5, v[j_act], x[j_act]);
} /* if (myRank == cur_rank) */ } /* if (myRank == cur_rank) */
} /* i */ } /* i */