diff --git a/TODO.md b/TODO.md index a81190a..2cceccc 100644 --- a/TODO.md +++ b/TODO.md @@ -1,7 +1,11 @@ TODO ==== -* Get rid of all global variables. +* Make the reamining arrays vectors or smart pointers. + +* Memory bug when reading HDF5? x and v not allocated. + +* Get rid of all global variables. Including config. * Break main() into smaller chunks; operations that are timed should become independent functions. diff --git a/black_holes.cpp b/black_holes.cpp index 4f5ffce..c14b970 100644 --- a/black_holes.cpp +++ b/black_holes.cpp @@ -99,7 +99,7 @@ void Black_hole_physics::adjust_post_newtonian( jrk2 += jrk2_corr; } -void Black_hole_physics::write_bh_data(double time_cur, double m[], double3 x[], double3 v[], const std::vector& pot, double3 a[], double3 adot[], double dt[]) +void Black_hole_physics::write_bh_data(double time_cur, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, double3 a[], double3 adot[], double dt[]) { // 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 diff --git a/black_holes.h b/black_holes.h index 3d15efa..402a628 100644 --- a/black_holes.h +++ b/black_holes.h @@ -46,7 +46,7 @@ public: const double dt_bh, // pn_usage should be const double3& acc1, double3& acc2, double3& jrk1, double3& jrk2); - void write_bh_data(double time_cur, double m[], double3 x[], double3 v[], const std::vector& pot, double3 a[], double3 adot[], double dt[]); + void write_bh_data(double time_cur, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, double3 a[], double3 adot[], double dt[]); public: //TODO make private double m1, m2; int count; @@ -58,11 +58,9 @@ public: //TODO make private Bbh_gravity bbh_grav; }; -//void write_bh_nb_data(int nb, int smbh_count, double time_cur, int N, double m[], double3 x[], double3 v[]); - class Write_bh_nb_data { public: - Write_bh_nb_data(int nb, int smbh_count, int N, double *m, double3 *x, double3 *v) + Write_bh_nb_data(int nb, int smbh_count, int N, const std::vector &m, const std::vector &x, const std::vector &v) : nb(nb), smbh_count(smbh_count), N(N), m(m), x(x), v(v) { ind_sort.resize(N); @@ -76,8 +74,8 @@ public: void operator()(double time_cur); private: int nb, smbh_count, N; - double *m; - double3 *x, *v; + const std::vector &m; + const std::vector &x, &v; std::vector ind_sort; std::vector var_sort; FILE *out; @@ -85,7 +83,7 @@ private: class Binary_smbh_influence_sphere_output { public: - Binary_smbh_influence_sphere_output(double factor, int N, double *m, double3 *x, double3 *v, const std::vector& pot, double *dt) + Binary_smbh_influence_sphere_output(double factor, int N, const std::vector &m, const std::vector &x, const std::vector &v, const std::vector &pot, double *dt) : factor(factor), m(m), x(x), v(v), pot(pot), dt(dt) { inf_event.assign(N, 0); @@ -99,9 +97,10 @@ public: void operator()(const std::vector& ind_act, int n_act, double timesteps, double time_cur); private: double factor; - const std::vector& pot; - double *m, /**pot,*/ *dt; - double3 *x, *v; + const std::vector &pot; + const std::vector &m; + double /**pot,*/ *dt; + const std::vector &x, &v; std::vector inf_event; FILE *out; }; \ No newline at end of file diff --git a/external.h b/external.h index aedcab5..08d989f 100644 --- a/external.h +++ b/external.h @@ -1,9 +1,10 @@ #pragma once +#include #include "double3.h" class External_gravity { public: - void apply(const int n_act, const double3 x[], const double3 v[], double pot[], double3 a[], double3 adot[]) + void apply(const int n_act, const std::vector &x, const std::vector &v, double pot[], double3 a[], double3 adot[]) { for (int i=0; iset_coordinates(x[i], v[i]); diff --git a/io.cpp b/io.cpp index 3a1b086..6bc8e5f 100644 --- a/io.cpp +++ b/io.cpp @@ -24,7 +24,7 @@ bool is_hdf5(std::string file_name) return result; } -void ascii_read(const std::string file_name, int& step_num, int& N, double& t, double **m, double3 **x, double3 **v) +void ascii_read(const std::string file_name, int &step_num, int &N, double& t, std::vector &m, std::vector &x, std::vector &v) { char rest[512]; int result; @@ -45,19 +45,19 @@ void ascii_read(const std::string file_name, int& step_num, int& N, double& t, d result = sscanf(str.c_str(), "%lf%s", &t, rest); if (result!=1) throw std::runtime_error("Error parsing line 3: expected one real number"); - *m = new double[N]; - *x = new double3[N]; - *v = new double3[N]; + m.resize(N); + x.resize(N); + v.resize(N); int i = -1; while (std::getline(file, str)) { if (++i > N) throw std::runtime_error("Error parsing line " + std::to_string(i+4) + ": particle out of range"); - result = sscanf(str.c_str(), "%*s %lf %lf %lf %lf %lf %lf %lf%s", &(*m)[i], &(*x)[i][0], &(*x)[i][1], &(*x)[i][2], &(*v)[i][0], &(*v)[i][1], &(*v)[i][2], rest); + result = sscanf(str.c_str(), "%*s %lf %lf %lf %lf %lf %lf %lf%s", &m[i], &(x[i].x), &(x[i].y), &(x[i].z), &(v[i].x), &(v[i].y), &(v[i].z), rest); } file.close(); } -void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, int precision=10) +void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector &m, const std::vector &x, const std::vector &v, int precision=10) { auto file = std::ofstream(file_name); if (!file.is_open()) throw std::runtime_error("Cannot open file for output"); @@ -146,7 +146,7 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub #endif } -void h5_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, const std::vector& 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, std::vector &m, std::vector &x, std::vector &v, const std::vector &pot, const double3 *acc, const double3 *jrk, const int extra_mode=0, const bool use_double_precision=true) { #ifdef HAS_HDF5 hid_t file_id, group_id, attribute_id, dataspace_id; @@ -174,9 +174,9 @@ void h5_write(const std::string file_name, const int step_num, const int N, cons H5Sclose(dataspace_id); }; - write_dataset("MASS", 1, (double*)m); - write_dataset("X", 2, (double*)x); - write_dataset("V", 2, (double*)v); + write_dataset("MASS", 1, m.data()); + write_dataset("X", 2, (double*)x.data()); + write_dataset("V", 2, (double*)v.data()); bool write_pot = (extra_mode ) & 1; bool write_acc = (extra_mode >> 1) & 1; diff --git a/io.h b/io.h index 8e8d2ad..874ca27 100644 --- a/io.h +++ b/io.h @@ -5,12 +5,12 @@ bool is_hdf5(std::string file_name); // This function is implemented independently of the HDF5 library -void ascii_read(const std::string file_name, int& step_num, int& N, double& t, double **m, double3 **x, double3 **v); +void ascii_read(const std::string file_name, int &step_num, int &N, double &t, std::vector &m, std::vector &x, std::vector &v); -void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, int precision=10); +void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const std::vector &m, const std::vector &x, const std::vector &v, int precision=10); 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 -void h5_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, const std::vector& pot, const double3 *acc, const double3 *jrk, const int write_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, std::vector &m, std::vector &x, std::vector &v, const std::vector &pot, const double3 *acc, const 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 diff --git a/phigrape.cpp b/phigrape.cpp index f8cc761..0d5b4dd 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -28,7 +28,6 @@ #include "grapite.h" #endif -Config *config; //chrono::steady_clock::time_point walltime_start; namespace std::chrono { @@ -58,7 +57,7 @@ public: acc_loc.resize(N); jrk_loc.resize(N); } - void operator()(const double t, const int n_act, std::vector& ind_act, const double3 x_act[], const double3 v_act[], + void operator()(const double t, const int n_act, std::vector &ind_act, std::vector &x_act, std::vector &v_act, std::vector& pot, double3 acc[], double3 jrk[]) { g6_set_ti(clusterid, t); @@ -84,7 +83,7 @@ private: std::vector acc_loc, jrk_loc; }; -void calc_ext_grav(std::vector &external_gravity_components, int n, double3 *x, double3 *v, double *pot, double3 *acc, double3* jrk) +void calc_ext_grav(std::vector &external_gravity_components, int n, const std::vector &x, const std::vector &v, double *pot, double3 *acc, double3* jrk) // TODO should just be a class that has this pointer array as a member { std::fill(pot, pot+n, 0.); @@ -94,7 +93,7 @@ void calc_ext_grav(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, double m[], double3 x[], double3 v[], const std::vector& 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 &m, const std::vector &x, const std::vector &v, const std::vector &pot, double pot_ext[]) { double E_pot = 0; for (int i=0; igrapite_active_search) { + if (grapite_active_search_flag) { min_t_loc = grapite_get_minimum_time(); } else #endif @@ -188,7 +187,7 @@ public: void get_active_indices(const double min_t, const double t[], const double dt[], int ind_act[], int& n_act) { #ifdef ETICS - if (config->grapite_active_search) { + if (grapite_active_search_flag) { int n_act_loc; grapite_active_search(min_t, ind_act_loc, &n_act_loc); if (myRank > 0) @@ -215,6 +214,7 @@ public: private: int myRank, n_proc, n_loc, N; int *ind_act_loc; + 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) @@ -254,6 +254,9 @@ inline double aarseth_step(const double eta, const double dt, const double3 a, c int main(int argc, char *argv[]) { + /* read the input parameters */ + const Config config("phigrape.conf"); + timer.start(); double timesteps=0.0, n_act_sum=0.0; @@ -280,23 +283,21 @@ int main(int argc, char *argv[]) /* Print the Rank and the names of processors */ printf("Rank of the processor %03d on %s \n", myRank, processor_name); - /* read the input parameters */ - config = new Config("phigrape.conf"); - int diskstep, N; double time_cur; // The memory for m, x, and v is allocated inside h5_read or ascii_read - double *m; - double3 *x, *v; - if (is_hdf5(config->input_file_name)) { + //double *m; + std::vector m; + std::vector x, v; + if (is_hdf5(config.input_file_name)) { #ifndef HAS_HDF5 fprintf(stderr, "ERROR: input file is in HDF5 format, but the code was compiled without HDF5 support\n"); return -1; #endif - h5_read(config->input_file_name, &diskstep, &N, &time_cur, m, x, v); + h5_read(config.input_file_name, &diskstep, &N, &time_cur, m.data(), x.data(), v.data()); } 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 ind(N); std::iota(begin(ind), end(ind), 0); @@ -308,16 +309,16 @@ int main(int argc, char *argv[]) int n_act; std::vector ind_act(N); std::vector pot_act_new(N); + std::vector x_act_new(N), v_act_new(N); double *pot_act_ext = new double[N]; - double3 *x_act_new = new double3[N], *v_act_new = new double3[N], - *a_act_new = new double3[N], *adot_act_new = new double3[N]; + double3 *a_act_new = new double3[N], *adot_act_new = new double3[N]; - double eps = config->eps; - double eta = config->eta; - double t_end = config->t_end; - double dt_disk = config->dt_disk; - double dt_contr = config->dt_contr; - double dt_bh = config->dt_bh; + double eps = config.eps; + double eta = config.eta; + double t_end = config.t_end; + double dt_disk = config.dt_disk; + double dt_contr = config.dt_contr; + double dt_bh = config.dt_bh; if (myRank == rootRank) { printf("\n"); @@ -333,11 +334,11 @@ int main(int argc, char *argv[]) if ((diskstep == 0) && (time_cur == 0)) { FILE *out = fopen("contr.dat", "w"); fclose(out); - if (config->live_smbh_output && (config->live_smbh_count > 0)) { + if (config.live_smbh_output && (config.live_smbh_count > 0)) { out = fopen("bh.dat", "w"); fclose(out); } - if ((config->live_smbh_neighbor_output) && (config->live_smbh_count > 0)) { + if ((config.live_smbh_neighbor_output) && (config.live_smbh_count > 0)) { out = fopen("bh_neighbors.dat", "w"); fclose(out); } @@ -345,16 +346,16 @@ int main(int argc, char *argv[]) } /* if (myRank == rootRank) */ double normalization_mass=1, normalization_length=1, normalization_velocity=1; - if (config->ext_units_physical) { - normalization_mass = 1/config->unit_mass; - normalization_length = 1000/config->unit_length; - normalization_velocity = 1.52484071426404437233e+01*sqrt(config->unit_length/config->unit_mass); + if (config.ext_units_physical) { + normalization_mass = 1/config.unit_mass; + normalization_length = 1000/config.unit_length; + normalization_velocity = 1.52484071426404437233e+01*sqrt(config.unit_length/config.unit_mass); } - Plummer ext_bulge(config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length); - Miyamoto_Nagai ext_disk(config->ext_m_disk*normalization_mass, config->ext_a_disk*normalization_length, config->ext_b_disk*normalization_length); - Plummer ext_halo_plummer(config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length); - Logarithmic_halo ext_log_halo(config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length); - Dehnen ext_dehnen(config->ext_dehnen_m*normalization_mass, config->ext_dehnen_r*normalization_length, config->ext_dehnen_gamma); + Plummer ext_bulge(config.ext_m_bulge*normalization_mass, config.ext_b_bulge*normalization_length); + Miyamoto_Nagai ext_disk(config.ext_m_disk*normalization_mass, config.ext_a_disk*normalization_length, config.ext_b_disk*normalization_length); + Plummer ext_halo_plummer(config.ext_m_halo_plummer*normalization_mass, config.ext_b_halo_plummer*normalization_length); + Logarithmic_halo ext_log_halo(config.ext_log_halo_v*normalization_velocity, config.ext_log_halo_r*normalization_length); + Dehnen ext_dehnen(config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma); std::vector external_gravity_components; external_gravity_components.push_back(&ext_bulge); @@ -371,11 +372,11 @@ int main(int argc, char *argv[]) } if (has_external_gravity) printf("External Potential: \n\n"); - if (ext_bulge.is_active) printf("m_bulge = %.4E b_bulge = %.4E\n", config->ext_m_bulge*normalization_mass, config->ext_b_bulge*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_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_bulge.is_active) printf("m_bulge = %.4E b_bulge = %.4E\n", config.ext_m_bulge*normalization_mass, config.ext_b_bulge*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_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); printf("\n"); fflush(stdout); @@ -399,13 +400,13 @@ int main(int argc, char *argv[]) /* some local settings for G6a boards */ int clusterid, numGPU; - if (config->devices_per_node==0) { + if (config.devices_per_node==0) { MPI_Comm shmcomm; MPI_Comm_split_type(MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, &shmcomm); MPI_Comm_size(shmcomm, &numGPU); MPI_Comm_rank(shmcomm, &clusterid); } else { - numGPU = config->devices_per_node; + numGPU = config.devices_per_node; clusterid = myRank % numGPU; } printf("Rank of the processor %03d : Number of GPUs %01d : Cluster ID %01d \n", myRank, numGPU, clusterid); @@ -418,30 +419,33 @@ int main(int argc, char *argv[]) g6_set_xunit(51); #ifdef ETICS - grapite_set_dev_exec_threshold(config->grapite_dev_exec_threshold); + grapite_set_dev_exec_threshold(config.grapite_dev_exec_threshold); + bool grapite_active_search_flag = config.grapite_active_search; +#else + bool grapite_active_search_flag = false; #endif 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); + 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) + if (config.live_smbh_count == 1) black_hole_physics = Black_hole_physics(m[0], 0, myRank, rootRank); - else if (config->live_smbh_count == 2) + 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()); + 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); + 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); + 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); + Write_bh_nb_data write_bh_nb_data(config.live_smbh_neighbor_number, config.live_smbh_count, N, m, x, v); #ifdef ETICS - grapite_read_particle_tags(N, config->grapite_mask_file_name.c_str(), myRank, n_loc); - grapite_set_dt_exp(config->dt_scf); + grapite_read_particle_tags(N, config.grapite_mask_file_name.c_str(), myRank, n_loc); + grapite_set_dt_exp(config.dt_scf); grapite_set_t_exp(time_cur); #endif @@ -470,10 +474,10 @@ int main(int argc, char *argv[]) /* define the all particles as a active on all the processors for the first time grav calc. */ calc_self_grav(time_cur, N, ind, x, v, pot, a, adot); - if (config->live_smbh_count == 2) { + if (config.live_smbh_count == 2) { 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->binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[0], 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]); } calc_ext_grav(external_gravity_components, N, x, v, pot_ext, a, adot); @@ -482,7 +486,7 @@ int main(int argc, char *argv[]) if (myRank == rootRank) energy_contr(time_cur, timesteps, n_act_sum, calc_self_grav.g6_calls, N, m, x, v, pot, pot_ext); #ifdef ETICS - if (config->etics_dump_coeffs && (diskstep==0)) { + if (config.etics_dump_coeffs && (diskstep==0)) { char out_fname[256]; sprintf(out_fname, "coeffs.%06d.%02d.dat", 0, myRank); grapite_dump(out_fname, 2); @@ -515,7 +519,7 @@ int main(int argc, char *argv[]) dt[j] = dt_tmp; - if (config->dt_min_warning && (myRank == 0)) { + if (config.dt_min_warning && (myRank == 0)) { if (dt[j] == dt_min) { printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[j], ind[j]); fflush(stdout); @@ -523,9 +527,9 @@ int main(int argc, char *argv[]) } } /* j */ - if (config->live_smbh_count > 0) { + if (config.live_smbh_count > 0) { double min_dt = *std::min_element(dt, dt+N); - for (int i=0; ilive_smbh_count; i++) dt[i] = min_dt; + for (int i=0; ilive_smbh_output) black_hole_physics.write_bh_data(time_cur, m, x, v, pot, a, adot, dt); + 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 (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur); } /* if (myRank == rootRank) */ @@ -558,8 +562,8 @@ int main(int argc, char *argv[]) /* Find the BH(s) indices in the active list */ int i_bh1=0, i_bh2=1; #ifdef ETICS - int n_bh = config->live_smbh_count; - if (config->grapite_active_search && (n_bh>0)) { + int n_bh = config.live_smbh_count; + if (config.grapite_active_search && (n_bh>0)) { int act_def_grapite_bh_count = 0; int i_bh[n_bh]; for (int i=0; ilive_smbh_count == 2) { + if (config.live_smbh_count == 2) { black_hole_physics.set_xv(x_act_new[i_bh1], x_act_new[i_bh2], v_act_new[i_bh1], v_act_new[i_bh2]); - if (config->live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(pot_act_new[i_bh1], pot_act_new[i_bh2], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]); - if (config->binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[i_bh1], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]); + if (config.live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(pot_act_new[i_bh1], pot_act_new[i_bh2], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]); + if (config.binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[i_bh1], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]); } /* Calculate gravity on active particles due to external forces */ @@ -608,7 +612,7 @@ int main(int argc, char *argv[]) //TODO make beautiful double eta_curr; - if ((config->live_smbh_count > 0) && (ind_act[i] < config->live_smbh_count)) eta_curr = eta_bh; + if ((config.live_smbh_count > 0) && (ind_act[i] < config.live_smbh_count)) eta_curr = eta_bh; else eta_curr = eta; double dt_new = aarseth_step(eta_curr, dt_tmp, a_act_new[i], adot_act_new[i], a2, a3); @@ -624,7 +628,7 @@ int main(int argc, char *argv[]) dt_tmp *= 2; } - if (config->dt_min_warning && (myRank == 0)) { + if (config.dt_min_warning && (myRank == 0)) { if (dt_tmp == dt_min) { printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt_tmp, ind_act[i]); fflush(stdout); @@ -643,12 +647,12 @@ int main(int argc, char *argv[]) } /* i */ /* define the min. dt over all the act. part. and set it also for the BH... */ - if (config->live_smbh_count > 0) { - if (config->live_smbh_count>=1) dt[0] = min_dt; - if (config->live_smbh_count==2) dt[1] = min_dt; + if (config.live_smbh_count > 0) { + if (config.live_smbh_count>=1) dt[0] = 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); } @@ -675,10 +679,10 @@ int main(int argc, char *argv[]) if (time_cur >= t_bh) { 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); + 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 (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur); } /* if (myRank == rootRank) */ @@ -689,7 +693,7 @@ int main(int argc, char *argv[]) if (myRank == rootRank) { energy_contr(time_cur, timesteps, n_act_sum, calc_self_grav.g6_calls, N, m, x, v, pot, pot_ext); /* write cont data */ - if (config->output_hdf5) h5_write("data.con", diskstep, N, time_cur, m, x, v, pot, a, adot, 0, true); + if (config.output_hdf5) h5_write("data.con", diskstep, N, time_cur, m, x, v, pot, a, adot, 0, true); else ascii_write("data.con", diskstep, N, time_cur, m, x, v, 16); } /* if (myRank == rootRank) */ @@ -717,12 +721,12 @@ int main(int argc, char *argv[]) diskstep++; if (myRank == rootRank) { sprintf(out_fname, "%06d", diskstep); - if (config->output_hdf5) h5_write(std::string(out_fname) + ".h5", diskstep, N, time_cur, m, x, v, pot, a, adot, config->output_extra_mode, config->output_hdf5_double_precision); - else ascii_write(std::string(out_fname) + ".dat", diskstep, N, time_cur, m, x, v, config->output_ascii_precision); + if (config.output_hdf5) h5_write(std::string(out_fname) + ".h5", diskstep, N, time_cur, m, x, v, pot, a, adot, config.output_extra_mode, config.output_hdf5_double_precision); + else ascii_write(std::string(out_fname) + ".dat", diskstep, N, time_cur, m, x, v, config.output_ascii_precision); } /* if (myRank == rootRank) */ #ifdef ETICS - if (config->etics_dump_coeffs) { + if (config.etics_dump_coeffs) { sprintf(out_fname, "coeffs.%06d.%02d.dat", diskstep, myRank); grapite_dump(out_fname, 2); } @@ -746,8 +750,7 @@ int main(int argc, char *argv[]) fflush(stdout); } /* if (myRank == rootRank) */ - delete config; - delete[] m; delete[] x; delete[] v; delete[] a; delete[] adot; delete[] pot_ext; delete[] t; delete[] dt; delete[] x_act_new; delete[] v_act_new; delete[] a_act_new; delete[] adot_act_new; delete[] pot_act_ext; + 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 */ MPI_Finalize();