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

Got rid of remaining manual allocations, moved some variable declarations to where they are needed

Browse source
This commit is contained in:
Yohai Meiron 2020-10-28 19:18:37 -04:00
parent aef0f26878
commit d9ed8e5131
8 changed files with 67 additions and 90 deletions
Download patch file Download diff file Expand all files Collapse all files

6
TODO.md
View file

@ -1,8 +1,10 @@
TODO TODO
==== ====
* Make the reamining arrays vectors or smart pointers.
* Memory bug when reading HDF5? x and v not allocated. * Memory bug when reading HDF5? x and v not allocated.
* Break main() into smaller chunks; operations that are timed should become independent functions. * Break main() into smaller chunks; operations that are timed should become independent functions.
* Const everything
* Remove unneeded includes

2
black_holes.cpp
View file

@ -99,7 +99,7 @@ void Black_hole_physics::adjust_post_newtonian(
jrk2 += jrk2_corr; jrk2 += jrk2_corr;
} }
void Black_hole_physics::write_bh_data(double time_cur, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, double3 a[], double3 adot[], double dt[]) void Black_hole_physics::write_bh_data(double time_cur, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &a, const std::vector<double3> &adot, const std::vector<double> &dt)
{ {
// This function logs data on the black hole(s). It uses both external data // This function logs data on the black hole(s). It uses both external data
// (the arguments to this function) and optionall internal data to this // (the arguments to this function) and optionall internal data to this

8
black_holes.h
View file

@ -46,7 +46,7 @@ public:
const double dt_bh, // pn_usage should be const const double dt_bh, // pn_usage should be const
double3& acc1, double3& acc2, double3& acc1, double3& acc2,
double3& jrk1, double3& jrk2); double3& jrk1, double3& jrk2);
void write_bh_data(double time_cur, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, double3 a[], double3 adot[], double dt[]); void write_bh_data(double time_cur, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &a, const std::vector<double3> &adot, const std::vector<double> &dt);
public: //TODO make private public: //TODO make private
double m1, m2; double m1, m2;
int count; int count;
@ -83,7 +83,7 @@ private:
class Binary_smbh_influence_sphere_output { class Binary_smbh_influence_sphere_output {
public: public:
Binary_smbh_influence_sphere_output(double factor, int N, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, double *dt) Binary_smbh_influence_sphere_output(double factor, int N, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double> &dt)
: factor(factor), m(m), x(x), v(v), pot(pot), dt(dt) : factor(factor), m(m), x(x), v(v), pot(pot), dt(dt)
{ {
inf_event.assign(N, 0); inf_event.assign(N, 0);
@ -97,9 +97,7 @@ public:
void operator()(const std::vector<int>& ind_act, int n_act, double timesteps, double time_cur); void operator()(const std::vector<int>& ind_act, int n_act, double timesteps, double time_cur);
private: private:
double factor; double factor;
const std::vector<double> &pot; const std::vector<double> &pot, &m, &dt;
const std::vector<double> &m;
double /**pot,*/ *dt;
const std::vector<double3> &x, &v; const std::vector<double3> &x, &v;
std::vector<int> inf_event; std::vector<int> inf_event;
FILE *out; FILE *out;

1
external.cpp
View file

@ -59,7 +59,6 @@ void Logarithmic_halo::calc_gravity()
jerk = tmp * (rv_ij * x - r2_r2_halo * v); jerk = tmp * (rv_ij * x - r2_r2_halo * v);
} }
void Dehnen::calc_gravity() void Dehnen::calc_gravity()
{ {

2
external.h
View file

@ -4,7 +4,7 @@
class External_gravity { class External_gravity {
public: public:
void apply(const int n_act, const std::vector<double3> &x, const std::vector<double3> &v, double pot[], double3 a[], double3 adot[]) void apply(const int n_act, const std::vector<double3> &x, const std::vector<double3> &v, std::vector<double> &pot, std::vector<double3> &a, std::vector<double3> &adot)
{ {
for (int i=0; i<n_act; i++) { for (int i=0; i<n_act; i++) {
this->set_coordinates(x[i], v[i]); this->set_coordinates(x[i], v[i]);

8
io.cpp
View file

@ -146,7 +146,7 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub
#endif #endif
} }
void h5_write(const std::string file_name, const int step_num, const int N, const double t, std::vector<double> &m, std::vector<double3> &x, std::vector<double3> &v, const std::vector<double> &pot, const double3 *acc, const double3 *jrk, const int extra_mode=0, const bool use_double_precision=true) void h5_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &acc, const std::vector<double3> &jrk, const int extra_mode=0, const bool use_double_precision=true)
{ {
#ifdef HAS_HDF5 #ifdef HAS_HDF5
hid_t file_id, group_id, attribute_id, dataspace_id; hid_t file_id, group_id, attribute_id, dataspace_id;
@ -174,7 +174,7 @@ void h5_write(const std::string file_name, const int step_num, const int N, cons
H5Sclose(dataspace_id); H5Sclose(dataspace_id);
}; };
write_dataset("MASS", 1, m.data()); write_dataset("MASS", 1, (double*)m.data()); // casting away const...
write_dataset("X", 2, (double*)x.data()); write_dataset("X", 2, (double*)x.data());
write_dataset("V", 2, (double*)v.data()); write_dataset("V", 2, (double*)v.data());
@ -182,8 +182,8 @@ void h5_write(const std::string file_name, const int step_num, const int N, cons
bool write_acc = (extra_mode >> 1) & 1; bool write_acc = (extra_mode >> 1) & 1;
bool write_jrk = (extra_mode >> 2) & 1; bool write_jrk = (extra_mode >> 2) & 1;
if (write_pot) write_dataset("POT", 1, (double*)pot.data()); if (write_pot) write_dataset("POT", 1, (double*)pot.data());
if (write_acc) write_dataset("ACC", 2, (double*)acc); if (write_acc) write_dataset("ACC", 2, (double*)acc.data());
if (write_jrk) write_dataset("JRK", 2, (double*)jrk); if (write_jrk) write_dataset("JRK", 2, (double*)jrk.data());
H5Gclose(group_id); H5Gclose(group_id);
H5Fclose(file_id); H5Fclose(file_id);

2
io.h
View file

@ -12,5 +12,5 @@ void ascii_write(const std::string file_name, const int step_num, const int N, c
void h5_read(const std::string file_name, int *step_num, int *N, double *t, double m[], double3 x[], double3 v[]); void h5_read(const std::string file_name, int *step_num, int *N, double *t, double m[], double3 x[], double3 v[]);
// In case the code is compiled without HDF5 support, the implementation of this function just throws an error // In case the code is compiled without HDF5 support, the implementation of this function just throws an error
void h5_write(const std::string file_name, const int step_num, const int N, const double t, std::vector<double> &m, std::vector<double3> &x, std::vector<double3> &v, const std::vector<double> &pot, const double3 *acc, const double3 *jrk, const int extra_mode=0, const bool use_double_precision=true); void h5_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double3> &acc, const std::vector<double3> &jrk, const int extra_mode=0, const bool use_double_precision=true);
// In case the code is compiled without HDF5 support, the implementation of this function just throws an error // In case the code is compiled without HDF5 support, the implementation of this function just throws an error

128
phigrape.cpp
View file

@ -58,21 +58,20 @@ public:
jrk_loc.resize(N); jrk_loc.resize(N);
} }
void operator()(const double t, const int n_act, std::vector<int> &ind_act, std::vector<double3> &x_act, std::vector<double3> &v_act, void operator()(const double t, const int n_act, std::vector<int> &ind_act, std::vector<double3> &x_act, std::vector<double3> &v_act,
std::vector<double>& pot, double3 acc[], double3 jrk[]) std::vector<double>& pot, std::vector<double3> &acc, std::vector<double3> &jrk)
{ {
g6_set_ti(clusterid, t); g6_set_ti(clusterid, t);
for (int i=0; i<n_act; i+=npipe) { for (int i=0; i<n_act; i+=npipe) {
int nn = npipe; int nn = npipe;
if (n_act-i < npipe) nn = n_act - i; if (n_act-i < npipe) nn = n_act - i;
//TODO any way we can clean up this ugly casting?
g6calc_firsthalf(clusterid, n_loc, nn, ind_act.data()+i, (double(*)[3])&x_act[i], (double(*)[3])&v_act[i], (double(*)[3])&acc_loc[i], (double(*)[3])&jrk_loc[i], &pot_loc[i], eps2, h2.data()); g6calc_firsthalf(clusterid, n_loc, nn, ind_act.data()+i, (double(*)[3])&x_act[i], (double(*)[3])&v_act[i], (double(*)[3])&acc_loc[i], (double(*)[3])&jrk_loc[i], &pot_loc[i], eps2, h2.data());
g6calc_lasthalf( clusterid, n_loc, nn, ind_act.data()+i, (double(*)[3])&x_act[i], (double(*)[3])&v_act[i], eps2, h2.data(), (double(*)[3])&acc_loc[i], (double(*)[3])&jrk_loc[i], &pot_loc[i]); g6calc_lasthalf( clusterid, n_loc, nn, ind_act.data()+i, (double(*)[3])&x_act[i], (double(*)[3])&v_act[i], eps2, h2.data(), (double(*)[3])&acc_loc[i], (double(*)[3])&jrk_loc[i], &pot_loc[i]);
g6_calls++; g6_calls++;
} /* i */ } /* i */
/* Reduce the "global" vectors from "local" on all the nodes */ /* Reduce the "global" vectors from "local" on all the nodes */
MPI_Allreduce(pot_loc.data(), pot.data(), n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(pot_loc.data(), pot.data(), n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(acc_loc.data(), acc, 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(acc_loc.data(), acc.data(), 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(jrk_loc.data(), jrk, 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(jrk_loc.data(), jrk.data(), 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
} }
double g6_calls; double g6_calls;
private: private:
@ -95,7 +94,7 @@ public:
if (component->is_active) return true; if (component->is_active) return true;
return false; return false;
} }
void operator()(int n, const std::vector<double3> &x, const std::vector<double3> &v, double *pot, double3 *acc, double3* jrk) void operator()(int n, const std::vector<double3> &x, const std::vector<double3> &v, std::vector<double> &pot, std::vector<double3> &acc, std::vector<double3> &jrk)
{ {
for (auto component : components) { for (auto component : components) {
if (component->is_active) if (component->is_active)
@ -106,7 +105,7 @@ private:
std::vector<External_gravity*> components; std::vector<External_gravity*> components;
}; };
void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, int N, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, double pot_ext[]) void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, int N, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, const std::vector<double> &pot_ext)
{ {
double E_pot = 0; double E_pot = 0;
for (int i=0; i<N; i++) E_pot += m[i]*pot[i]; for (int i=0; i<N; i++) E_pot += m[i]*pot[i];
@ -175,10 +174,9 @@ public:
Active_search(const int myRank, const int n_proc, const int n_loc, const int N, bool grapite_active_search_flag) Active_search(const int myRank, const int n_proc, const int n_loc, const int N, bool grapite_active_search_flag)
: myRank(myRank), n_proc(n_proc), n_loc(n_loc), N(N), grapite_active_search_flag(grapite_active_search_flag) : myRank(myRank), n_proc(n_proc), n_loc(n_loc), N(N), grapite_active_search_flag(grapite_active_search_flag)
{ {
ind_act_loc = new int[n_loc]; ind_act_loc.resize(n_loc);
} }
~Active_search() { delete[] ind_act_loc; }; double get_minimum_time(const std::vector<double> &t, std::vector<double> &dt)
double get_minimum_time(const double t[], const double dt[])
{ {
double min_t_loc, min_t; double min_t_loc, min_t;
#ifdef ETICS #ifdef ETICS
@ -197,7 +195,7 @@ public:
MPI_Allreduce(&min_t_loc, &min_t, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); MPI_Allreduce(&min_t_loc, &min_t, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
return min_t; return min_t;
} }
void get_active_indices(const double min_t, const double t[], const double dt[], int ind_act[], int& n_act) void get_active_indices(const double min_t, const std::vector<double> &t, const std::vector<double> &dt, std::vector<int> &ind_act, int &n_act)
{ {
#ifdef ETICS #ifdef ETICS
if (grapite_active_search_flag) { if (grapite_active_search_flag) {
@ -226,11 +224,11 @@ public:
} }
private: private:
int myRank, n_proc, n_loc, N; int myRank, n_proc, n_loc, N;
int *ind_act_loc; std::vector<int> ind_act_loc;
bool grapite_active_search_flag; bool grapite_active_search_flag;
}; };
inline void calc_high_derivatives(const double dt_tmp, const double3 a_old, const double3 a_new, const double3 a1_old, const double3 a1_new, double3& a2, double3& a3) inline void calc_high_derivatives(const double dt_tmp, const double3 &a_old, const double3 &a_new, const double3 &a1_old, const double3 &a1_new, double3 &a2, double3 &a3)
{ {
double dtinv = 1/dt_tmp; double dtinv = 1/dt_tmp;
double dt2inv = dtinv*dtinv; double dt2inv = dtinv*dtinv;
@ -244,7 +242,7 @@ inline void calc_high_derivatives(const double dt_tmp, const double3 a_old, cons
a3 = 12*a0mia1*dt3inv + 6*ad0plad1*dt2inv; a3 = 12*a0mia1*dt3inv + 6*ad0plad1*dt2inv;
} }
inline void corrector(const double dt_tmp, const double3 a2, const double3 a3, double3& x, double3& v) inline void corrector(const double dt_tmp, const double3 &a2, const double3 &a3, double3 &x, double3 &v)
{ {
double dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0; double dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0;
double dt4over24 = dt3over6*dt_tmp/4.0; double dt4over24 = dt3over6*dt_tmp/4.0;
@ -254,7 +252,7 @@ inline void corrector(const double dt_tmp, const double3 a2, const double3 a3, d
v += dt3over6*a2 + dt4over24*a3; v += dt3over6*a2 + dt4over24*a3;
} }
inline double aarseth_step(const double eta, const double dt, const double3 a, const double3 a1, const double3 a2, const double3 a3) inline double aarseth_step(const double eta, const double dt, const double3 &a, const double3 &a1, const double3 &a2, const double3 &a3)
{ {
double a1abs = a.norm(); double a1abs = a.norm();
double adot1abs = a1.norm(); double adot1abs = a1.norm();
@ -274,8 +272,6 @@ int main(int argc, char *argv[])
double timesteps=0.0, n_act_sum=0.0; double timesteps=0.0, n_act_sum=0.0;
double3 zeros = {0, 0, 0}; // Dummy; can't really be const because of the GRAPE interface.
/* INIT the rand() !!! */ /* INIT the rand() !!! */
srand(19640916); /* it is just my birthday :-) */ srand(19640916); /* it is just my birthday :-) */
@ -298,8 +294,6 @@ int main(int argc, char *argv[])
int diskstep, N; int diskstep, N;
double time_cur; double time_cur;
// The memory for m, x, and v is allocated inside h5_read or ascii_read
//double *m;
std::vector<double> m; std::vector<double> m;
std::vector<double3> x, v; std::vector<double3> x, v;
if (is_hdf5(config.input_file_name)) { if (is_hdf5(config.input_file_name)) {
@ -312,20 +306,6 @@ int main(int argc, char *argv[])
else else
ascii_read(config.input_file_name, diskstep, N, time_cur, m, x, v); ascii_read(config.input_file_name, diskstep, N, time_cur, m, x, v);
std::vector<int> ind(N);
std::iota(begin(ind), end(ind), 0);
double3 *a = new double3[N], *adot = new double3[N];
std::vector<double> pot(N);
double *pot_ext = new double[N], *t = new double[N], *dt = new double[N];
/* data for active particles */
int n_act;
std::vector<int> ind_act(N);
std::vector<double> pot_act_new(N);
std::vector<double3> x_act_new(N), v_act_new(N);
double *pot_act_ext = new double[N];
double3 *a_act_new = new double3[N], *adot_act_new = new double3[N];
double eps = config.eps; double eps = config.eps;
double eta = config.eta; double eta = config.eta;
double t_end = config.t_end; double t_end = config.t_end;
@ -382,8 +362,7 @@ int main(int argc, char *argv[])
if (ext_disk.is_active) printf("m_disk = %.4E a_disk = %.4E b_disk = %.4E\n", config.ext_m_disk*normalization_mass, config.ext_a_disk*normalization_length, config.ext_b_disk*normalization_length); if (ext_disk.is_active) printf("m_disk = %.4E a_disk = %.4E b_disk = %.4E\n", config.ext_m_disk*normalization_mass, config.ext_a_disk*normalization_length, config.ext_b_disk*normalization_length);
if (ext_halo_plummer.is_active) printf("m_halo = %.4E b_halo = %.4E\n", config.ext_m_halo_plummer*normalization_mass, config.ext_b_halo_plummer*normalization_length); if (ext_halo_plummer.is_active) printf("m_halo = %.4E b_halo = %.4E\n", config.ext_m_halo_plummer*normalization_mass, config.ext_b_halo_plummer*normalization_length);
if (ext_log_halo.is_active) printf("v_halo = %.6E r_halo = %.6E \n", config.ext_log_halo_v*normalization_velocity, config.ext_log_halo_r*normalization_length); if (ext_log_halo.is_active) printf("v_halo = %.6E r_halo = %.6E \n", config.ext_log_halo_v*normalization_velocity, config.ext_log_halo_r*normalization_length);
if (ext_dehnen.is_active) printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n", config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma); if (ext_dehnen.is_active) printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n\n", config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma);
printf("\n");
fflush(stdout); fflush(stdout);
double eta_bh = eta/ETA_BH_CORR; double eta_bh = eta/ETA_BH_CORR;
@ -401,9 +380,6 @@ int main(int argc, char *argv[])
fflush(stdout); fflush(stdout);
} /* if (myRank == rootRank) */ } /* if (myRank == rootRank) */
std::fill(t, t+N, time_cur);
std::fill(dt, dt+N, dt_min);
/* some local settings for G6a boards */ /* some local settings for G6a boards */
int clusterid, numGPU; int clusterid, numGPU;
if (config.devices_per_node==0) { if (config.devices_per_node==0) {
@ -431,33 +407,19 @@ int main(int argc, char *argv[])
#endif #endif
int n_loc = N/n_proc; int n_loc = N/n_proc;
Calc_self_grav calc_self_grav(N, n_loc, clusterid, npipe, eps);
Active_search active_search(myRank, n_proc, n_loc, N, grapite_active_search_flag);
Black_hole_physics black_hole_physics;
if (config.live_smbh_count == 1)
black_hole_physics = Black_hole_physics(m[0], 0, myRank, rootRank);
else if (config.live_smbh_count == 2)
black_hole_physics = Black_hole_physics(m[0], m[1], myRank, rootRank);
if (config.binary_smbh_pn) {
black_hole_physics.set_post_newtonian(config.pn_c, config.pn_usage.data());
if (config.pn_usage[6]) black_hole_physics.set_spins(config.smbh1_spin.data(), config.smbh2_spin.data());
}
black_hole_physics.set_softening(eps, config.live_smbh_custom_eps);
Binary_smbh_influence_sphere_output binary_smbh_influence_sphere_output(config.binary_smbh_influence_radius_factor, N, m, x, v, pot, dt);
Write_bh_nb_data write_bh_nb_data(config.live_smbh_neighbor_number, config.live_smbh_count, N, m, x, v);
#ifdef ETICS #ifdef ETICS
grapite_read_particle_tags(N, config.grapite_mask_file_name.c_str(), myRank, n_loc); grapite_read_particle_tags(N, config.grapite_mask_file_name.c_str(), myRank, n_loc);
grapite_set_dt_exp(config.dt_scf); grapite_set_dt_exp(config.dt_scf);
grapite_set_t_exp(time_cur); grapite_set_t_exp(time_cur);
#endif #endif
std::vector<int> ind(N);
std::iota(begin(ind), end(ind), 0);
/* load the nj particles to the G6 */ /* load the nj particles to the G6 */
double3 zeros = {0, 0, 0}; // Dummy; can't really be const because of the GRAPE interface.
for (int k=0; k<n_loc; k++) { for (int k=0; k<n_loc; k++) {
int j = k + myRank*n_loc; int j = k + myRank*n_loc;
g6_set_j_particle(clusterid, k, ind[j], t[j], dt[j], m[j], zeros, zeros, zeros, v[j], x[j]); g6_set_j_particle(clusterid, k, ind[j], time_cur, dt_min, m[j], zeros, zeros, zeros, v[j], x[j]);
} /* k */ } /* k */
#ifdef ETICS #ifdef ETICS
@ -476,16 +438,29 @@ int main(int argc, char *argv[])
} }
#endif #endif
std::vector<double3> a(N), adot(N);
std::vector<double> pot(N);
/* define the all particles as a active on all the processors for the first time grav calc. */ /* define the all particles as a active on all the processors for the first time grav calc. */
Calc_self_grav calc_self_grav(N, n_loc, clusterid, npipe, eps);
calc_self_grav(time_cur, N, ind, x, v, pot, a, adot); calc_self_grav(time_cur, N, ind, x, v, pot, a, adot);
if (config.live_smbh_count == 2) { Black_hole_physics black_hole_physics;
if (config.live_smbh_count == 1)
black_hole_physics = Black_hole_physics(m[0], 0, myRank, rootRank);
else if (config.live_smbh_count == 2) {
black_hole_physics = Black_hole_physics(m[0], m[1], myRank, rootRank);
black_hole_physics.set_xv(x[0], x[1], v[0], v[1]); black_hole_physics.set_xv(x[0], x[1], v[0], v[1]);
if (config.live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(pot[0], pot[1], a[0], a[1], adot[0], adot[1]); if (config.live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(pot[0], pot[1], a[0], a[1], adot[0], adot[1]);
if (config.binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[0], a[0], a[1], adot[0], adot[1]); if (config.binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt_min, a[0], a[1], adot[0], adot[1]);
} }
if (config.binary_smbh_pn) {
black_hole_physics.set_post_newtonian(config.pn_c, config.pn_usage.data());
if (config.pn_usage[6]) black_hole_physics.set_spins(config.smbh1_spin.data(), config.smbh2_spin.data());
}
black_hole_physics.set_softening(eps, config.live_smbh_custom_eps);
std::fill(pot_ext, pot_ext+N, 0.); std::vector<double> pot_ext(N, 0.);
calc_ext_grav(N, x, v, pot_ext, a, adot); calc_ext_grav(N, x, v, pot_ext, a, adot);
/* Energy control... */ /* Energy control... */
@ -507,6 +482,7 @@ int main(int argc, char *argv[])
} }
#endif #endif
std::vector<double> dt(N);
/* Define initial timestep for all particles on all nodes */ /* Define initial timestep for all particles on all nodes */
for (int j=0; j<N; j++) { for (int j=0; j<N; j++) {
double a2_mod = a[j].norm2(); double a2_mod = a[j].norm2();
@ -534,28 +510,33 @@ int main(int argc, char *argv[])
} /* j */ } /* j */
if (config.live_smbh_count > 0) { if (config.live_smbh_count > 0) {
double min_dt = *std::min_element(dt, dt+N); double min_dt = *std::min_element(begin(dt), end(dt));
for (int i=0; i<config.live_smbh_count; i++) dt[i] = min_dt; for (int i=0; i<config.live_smbh_count; i++) dt[i] = min_dt;
} }
/* load the new values for particles to the local GRAPEs */ /* load the new values for particles to the local GRAPEs */
for (int k=0; k<n_loc; k++) { for (int k=0; k<n_loc; k++) {
int j = k + myRank*n_loc; int j = k + myRank*n_loc;
g6_set_j_particle(clusterid, k, ind[j], t[j], dt[j], m[j], zeros, adot[j]*(1./6.), a[j]*0.5, v[j], x[j]); g6_set_j_particle(clusterid, k, ind[j], time_cur, dt[j], m[j], zeros, adot[j]*(1./6.), a[j]*0.5, v[j], x[j]);
} /* k */ } /* k */
if (myRank == rootRank) {
/* Write BH data... */
if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, m, x, v, pot, a, adot, dt);
/* Write BH NB data... */
if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur);
} /* if (myRank == rootRank) */
timesteps = 0.0; // Why won't those two be long long instead of double + should include the zeroth step timesteps = 0.0; // Why won't those two be long long instead of double + should include the zeroth step
n_act_sum = 0.0; n_act_sum = 0.0;
std::vector<int> ind_act(N);
std::vector<double3> x_act_new(N), v_act_new(N), a_act_new(N), adot_act_new(N);
std::vector<double> t(N, time_cur), pot_act_new(N);
// Functors for the main integration loop
Active_search active_search(myRank, n_proc, n_loc, N, grapite_active_search_flag);
Binary_smbh_influence_sphere_output binary_smbh_influence_sphere_output(config.binary_smbh_influence_radius_factor, N, m, x, v, pot, dt);
Write_bh_nb_data write_bh_nb_data(config.live_smbh_neighbor_number, config.live_smbh_count, N, m, x, v);
if (myRank == rootRank) {
if (config.live_smbh_output) black_hole_physics.write_bh_data(time_cur, m, x, v, pot, a, adot, dt);
if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur);
} /* if (myRank == rootRank) */
/* The main integration loop */ /* The main integration loop */
while (time_cur <= t_end) { while (time_cur <= t_end) {
@ -563,7 +544,8 @@ int main(int argc, char *argv[])
double min_t = active_search.get_minimum_time(t, dt); double min_t = active_search.get_minimum_time(t, dt);
/* Get indices of all particles that will be active in this bunch */ /* Get indices of all particles that will be active in this bunch */
active_search.get_active_indices(min_t, t, dt, ind_act.data(), n_act); int n_act;
active_search.get_active_indices(min_t, t, dt, ind_act, n_act);
/* 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;
@ -603,7 +585,7 @@ int main(int argc, char *argv[])
} }
/* Calculate gravity on active particles due to external forces */ /* Calculate gravity on active particles due to external forces */
std::fill(pot_act_ext, pot_act_ext+n_act, 0.); std::vector<double> pot_act_ext(N, 0.);
calc_ext_grav(n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new); calc_ext_grav(n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new);
/* correct the active particles positions etc... on all the nodes */ /* correct the active particles positions etc... on all the nodes */
@ -757,10 +739,6 @@ int main(int argc, char *argv[])
fflush(stdout); fflush(stdout);
} /* if (myRank == rootRank) */ } /* if (myRank == rootRank) */
delete[] a; delete[] adot; delete[] pot_ext; delete[] t; delete[] dt; delete[] a_act_new; delete[] adot_act_new; delete[] pot_act_ext;
/* Finalize the MPI work */ /* Finalize the MPI work */
MPI_Finalize(); MPI_Finalize();
return 0;
} }