diff --git a/config.cpp b/config.cpp index cce97cd..36d2b16 100644 --- a/config.cpp +++ b/config.cpp @@ -72,7 +72,6 @@ template<> bool Config::string_cast(const std::string str) throw std::runtime_error("Cannot convert \"" + str + "\" into a bool"); } -#include template<> std::vector Config::string_cast(const std::string str) { auto error = std::runtime_error("Cannot convert \"" + str + "\" into an integer array"); @@ -90,6 +89,23 @@ template<> std::vector Config::string_cast(const std::string str) return result; } +template<> std::vector Config::string_cast(const std::string str) +{ + auto error = std::runtime_error("Cannot convert \"" + str + "\" into an integer array"); + if (!( (str.front()=='{') && (str.back()=='}'))) throw error; + std::string new_str = strip(str.substr(1, str.length()-2)); + std::replace(new_str.begin(), new_str.end(), ',', ' '); + + std::vector result; + while (new_str.length() > 0) { + size_t idx; + auto value = std::stod(new_str, &idx); + result.push_back(value); + new_str = new_str.substr(idx, new_str.length()-idx); + } + return result; +} + // For mandatory parameters template @@ -109,6 +125,8 @@ T Config::get_parameter(Dictionary dictionary, std::string name, T default_value else return string_cast((*item).second); } +#include + void Config::error_checking() { #ifndef HAS_HDF5 @@ -130,9 +148,21 @@ void Config::error_checking() throw std::runtime_error("Binary black hole influence sphere output (binary_smbh_influence_sphere_output=true) requires exactly two live black holes (live_smbh_count=2)"); if (pn_usage.size() != 7) throw std::runtime_error("PN usage array (pn_usage) must have exactly seven components"); - for (int i=0; i<7; i++) - if (!((pn_usage[i] == 0) || (pn_usage[i] == 1))) - throw std::runtime_error("PN usage array (pn_usage) must have ones and zeros only"); + if (binary_smbh_pn) + for (int i=0; i<7; i++) + if (!((pn_usage[i] == 0) || (pn_usage[i] == 1))) + throw std::runtime_error("PN usage array (pn_usage) must be a 7-component vector filled with ones and zeros only"); + if ((smbh1_spin.size()!=3) || (smbh2_spin.size()!=3)) + throw std::runtime_error("Spins must be three-component vectors"); + if ((pn_usage[6]==1) && ((smbh1_spin[0]==nix) || (smbh2_spin[0]==nix))) + throw std::runtime_error("Please define smbh1_spin and smbh2_spin or disable the spin by setting the last component of pn_usage to zero"); + std::cout << smbh1_spin[0] << std::endl; + std::cout << smbh1_spin[1] << std::endl; + std::cout << smbh1_spin[2] << std::endl; + if ((pn_usage[6]==0) && ((smbh1_spin[0]!=nix) || (smbh2_spin[0]!=nix))) + throw std::runtime_error("Spins (smbh1_spin and smbh2_spin) may not be defined if the last element of pn_usage is set to zero"); + + if (ext_units_physical && ((unit_mass == 0) || (unit_length == 0))) throw std::runtime_error("Physical units for external gravity (ext_units_physical) requires ext_unit_mass and ext_unit_length to be positive numbers"); if ((ext_m_bulge > 0) && (ext_b_bulge < 0)) @@ -177,8 +207,10 @@ Config::Config(std::string file_name) binary_smbh_influence_radius_factor = get_parameter(dictionary, "binary_smbh_influence_radius_factor", 10.); binary_smbh_pn = get_parameter(dictionary, "binary_smbh_pn", false); - pn_usage = get_parameter>(dictionary, "pn_usage", std::vector({1,1,1,1,1,1,1})); + pn_usage = get_parameter>(dictionary, "pn_usage", std::vector({-1,-1,-1,-1,-1,-1,-1})); pn_c = get_parameter(dictionary, "pn_c", 0); + smbh1_spin = get_parameter>(dictionary, "smbh1_spin", std::vector({nix,nix,nix})); + smbh2_spin = get_parameter>(dictionary, "smbh2_spin", std::vector({nix,nix,nix})); ext_units_physical = get_parameter(dictionary, "ext_units_physical", false); unit_mass = get_parameter(dictionary, "unit_mass", !ext_units_physical); diff --git a/config.h b/config.h index f331315..6805646 100644 --- a/config.h +++ b/config.h @@ -7,48 +7,50 @@ class Config { public: Config(std::string file_name); - double eps; - double t_end; - double dt_disk; - double dt_contr; - double dt_bh; - double eta; - std::string input_file_name; - int devices_per_node; - bool dt_min_warning; + double eps; + double t_end; + double dt_disk; + double dt_contr; + double dt_bh; + double eta; + std::string input_file_name; + int devices_per_node; + bool dt_min_warning; + + bool output_hdf5; + bool output_hdf5_double_precision; + int output_ascii_precision; + int output_extra_mode; + + int live_smbh_count; + double live_smbh_custom_eps; + bool live_smbh_output; + bool live_smbh_neighbor_output; + int live_smbh_neighbor_number; + + bool binary_smbh_pn; + bool binary_smbh_influence_sphere_output; + double binary_smbh_influence_radius_factor; + std::vector pn_usage; + double pn_c; + std::vector smbh1_spin; + std::vector smbh2_spin; - bool output_hdf5; - bool output_hdf5_double_precision; - int output_ascii_precision; - int output_extra_mode; - - int live_smbh_count; - double live_smbh_custom_eps; - bool live_smbh_output; - bool live_smbh_neighbor_output; - int live_smbh_neighbor_number; - - bool binary_smbh_pn; - bool binary_smbh_influence_sphere_output; - double binary_smbh_influence_radius_factor; - std::vector pn_usage; - double pn_c; - - bool ext_units_physical; - double unit_mass; - double unit_length; - double ext_m_bulge; - double ext_b_bulge; - double ext_m_disk; - double ext_a_disk; - double ext_b_disk; - double ext_m_halo_plummer; - double ext_b_halo_plummer; - double ext_log_halo_r; - double ext_log_halo_v; - double ext_dehnen_m; - double ext_dehnen_r; - double ext_dehnen_gamma; + bool ext_units_physical; + double unit_mass; + double unit_length; + double ext_m_bulge; + double ext_b_bulge; + double ext_m_disk; + double ext_a_disk; + double ext_b_disk; + double ext_m_halo_plummer; + double ext_b_halo_plummer; + double ext_log_halo_r; + double ext_log_halo_v; + double ext_dehnen_m; + double ext_dehnen_r; + double ext_dehnen_gamma; private: using Dictionary = std::unordered_map; diff --git a/init.py b/init.py index a12d55c..649f75f 100644 --- a/init.py +++ b/init.py @@ -13,7 +13,7 @@ def gen_plum(N, seed=None, RMAX=10): X = np.sqrt(R**2 - Z**2) * np.cos(2*np.pi*X3) Y = np.sqrt(R**2 - Z**2) * np.sin(2*np.pi*X3) - Ve = np.sqrt(2)*(1.0 + R**2)**(-0.25); + Ve = np.sqrt(2)*(1.0 + R**2)**(-0.25) X4, X5 = 0, 0 while 0.1*X5 >= X4**2*(1-X4**2)**3.5: @@ -22,9 +22,9 @@ def gen_plum(N, seed=None, RMAX=10): V = Ve*X4 X6, X7 = np.random.random(2) - Vz = (1 - 2*X6)*V; - Vx = np.sqrt(V**2 - Vz**2) * np.cos(2*np.pi*X7); - Vy = np.sqrt(V**2 - Vz**2) * np.sin(2*np.pi*X7); + Vz = (1 - 2*X6)*V + Vx = np.sqrt(V**2 - Vz**2) * np.cos(2*np.pi*X7) + Vy = np.sqrt(V**2 - Vz**2) * np.sin(2*np.pi*X7) X, Y, Z = np.array([X, Y, Z])*3*np.pi/16 Vx, Vy, Vz = np.array([Vx, Vy, Vz])/np.sqrt(3*np.pi/16) @@ -38,7 +38,7 @@ def kepler_to_cartesian(a, e, i, Omega, w, nu, G=1.0, M=1.0): result = [] for arg in args: result.append(np.atleast_1d(arg)) return result - a, e, i, Omega, w, nu = to_arrays(a, e, i, Omega, w, nu) + a, e, i, Omega, w, nu = to_arrays(a, e, i, Omega, w, nu) # pylint: disable=unbalanced-tuple-unpacking P = [np.cos(w)*np.cos(Omega) - np.sin(w)*np.cos(i)*np.sin(Omega), np.cos(w)*np.sin(Omega) + np.sin(w)*np.cos(i)*np.cos(Omega), np.sin(w)*np.sin(i)] diff --git a/phigrape.cpp b/phigrape.cpp index f4b5267..eae49f2 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -594,80 +594,79 @@ inline double aarseth_step(const double eta, const double dt_tmp, const double3 return sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); } -void binary_smbh_influence_sphere_output(int i_bh1, int i_bh2, int ind_act[], double m_act[], double3 x_act[], double3 v_act[], double pot_act[], double dt_act[], int n_act, double timesteps, double time_cur, double factor, int inf_event[]) +void binary_smbh_influence_sphere_output(int i_bh1, int i_bh2, int ind_act[], int n_act, double timesteps, double time_cur, double factor, int inf_event[], Source_particle_list *source_particle_list) { //TODO !!IMPORTANT!! only open the file IF THERE IS ANYTHING TO WRITE!!! //TODO inf_event to be static or something? auto out = fopen("bbh_inf.dat", "a"); - double m_bh1 = m_act[i_bh1]; - double m_bh2 = m_act[i_bh2]; - - double3 x_bh1 = x_act[i_bh1]; - double3 x_bh2 = x_act[i_bh2]; - - double3 v_bh1 = v_act[i_bh1]; - double3 v_bh2 = v_act[i_bh2]; + double m_bh1 = source_particle_list->m[0]; + double m_bh2 = source_particle_list->m[1]; + double3 x_bh1 = source_particle_list->x[0]; + double3 x_bh2 = source_particle_list->x[1]; + double3 v_bh1 = source_particle_list->v[0]; + double3 v_bh2 = source_particle_list->v[1]; double3 x_bbhc = (m_bh1*x_bh1 + m_bh2*x_bh2)/(m_bh1 + m_bh2); double3 v_bbhc = (m_bh1*v_bh1 + m_bh2*v_bh2)/(m_bh1 + m_bh2); - double DR2 = SQR(x_bh1[0] - x_bh2[0]) + SQR(x_bh1[1] - x_bh2[1]) + SQR(x_bh1[2] - x_bh2[2]); - double DV2 = SQR(v_bh1[0] - v_bh2[0]) + SQR(v_bh1[1] - v_bh2[1]) + SQR(v_bh1[2] - v_bh2[2]); - + double DR2 = (x_bh1 - x_bh2).norm2(); + double DV2 = (v_bh1 - v_bh2).norm2(); double EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2; - double SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB; double SEMI_a2 = SQR(SEMI_a); - for (int i=2; ipot[0]; + double& pot_bh2 = source_particle_list->pot[1]; + double& m_act = source_particle_list->m[j_act]; + double3& x_act = source_particle_list->x[j_act]; + double3& v_act = source_particle_list->v[j_act]; + double& dt_act = source_particle_list->dt[j_act]; + double& pot_act = source_particle_list->pot[j_act]; + double tmp_r2 = (x_act - x_bbhc).norm2(); if (tmp_r2 < SEMI_a2*SQR(factor)) { - if (inf_event[iii] == 0) { + if (inf_event[j_act] == 0) { fprintf(out,"INF1 %.6E %.16E %07d %07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E \n", - timesteps, time_cur, i, ind_act[i], + timesteps, time_cur, i, j_act, sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2], - m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0], - m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1], + m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_bh1, + m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_bh2, sqrt(tmp_r2), - m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i], - dt_act[i]); + m_act, x_act[0], x_act[1], x_act[2], v_act[0], v_act[1], v_act[2], pot_act, + dt_act); - inf_event[iii] = 1; + inf_event[j_act] = 1; } } else { - if (inf_event[iii] == 1) { + if (inf_event[j_act] == 1) { fprintf(out,"INF2 %.6E %.16E %07d %07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E \n", timesteps, time_cur, i, ind_act[i], sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2], - m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0], - m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1], + m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_bh1, + m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_bh2, sqrt(tmp_r2), - m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i], - dt_act[i]); + m_act, x_act[0], x_act[1], x_act[2], v_act[0], v_act[1], v_act[2], pot_act, + dt_act); } - inf_event[iii] = 0; + inf_event[j_act] = 0; } /* if (tmp_r2 < DR2*R_INF2) */ } /* i */ fclose(out); - } -void adjust_bsmbh_softening(const double eps, const double eps_bh, const int i_bh1, const int i_bh2, const double m_act[], const double3 x_act_new[], const double3 v_act_new[], double pot_act_new[], double3 a_act_new[], double3 adot_act_new[]) +inline void adjust_bsmbh_softening(const double eps, const double eps_bh, const int i_bh1, const int i_bh2, const double m_bh1, const double m_bh2, const double3 x_act_new[], const double3 v_act_new[], double pot_act_new[], double3 a_act_new[], double3 adot_act_new[]) { - double m_bh1 = m_act[i_bh1]; - double m_bh2 = m_act[i_bh2]; - double3 x_bh1 = x_act_new[i_bh1]; double3 v_bh1 = v_act_new[i_bh1]; @@ -767,9 +766,6 @@ int main(int argc, char *argv[]) FILE *out; - double C_NB; - int usedOrNot[7]; - double pot_ext[N_MAX], pot_act_ext[N_MAX]; // for all EXTPOT int i_bh1, i_bh2; @@ -779,9 +775,6 @@ int main(int argc, char *argv[]) int inf_event[N_MAX]; double3 x_bbhc, v_bbhc; - double s_bh1[3] = {0.0, 0.0, 1.0}; - double s_bh2[3] = {0.0, 0.0, 1.0}; - double3 zeros = {0, 0, 0}; // Dummy; can't really be const because of the GRAPE interface. /* INIT the rand() !!! */ @@ -808,8 +801,6 @@ int main(int argc, char *argv[]) /* read the input parameters to the rootRank */ config = new Config("phigrape.conf"); - C_NB = config->pn_c; - std::copy(config->pn_usage.begin(), config->pn_usage.end(), usedOrNot); if (is_hdf5(config->input_file_name)) { #ifndef HAS_HDF5 @@ -1006,10 +997,10 @@ int main(int argc, char *argv[]) /* load the nj particles to the G6 */ - for (int j=0; j BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk - dt_bh_tmp = dt[0]; - - tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1, - m_bh2, x_bh2, v_bh2, s_bh2, - C_NB, dt_bh_tmp, usedOrNot, - (double(*)[3])a_pn1, (double(*)[3])adot_pn1, - (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank); + tmp_i = calc_force_pn_BH(m[i_bh1], x[i_bh1], v[i_bh1], config->smbh1_spin.data(), + m[i_bh2], x[i_bh2], v[i_bh2], config->smbh2_spin.data(), + config->pn_c, dt[i_bh1], config->pn_usage.data(), + (double(*)[3])a_pn1, (double(*)[3])adot_pn1, + (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank); a[i_bh1] += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6]; a[i_bh2] += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6]; @@ -1187,14 +1176,12 @@ int main(int argc, char *argv[]) /* Define initial timestep for all particles on all nodes */ - for (int i=0; i dt_max) dt_tmp = dt_max; - dt[i] = dt_tmp; + dt[j] = dt_tmp; if (config->dt_min_warning && (myRank == 0)) { - if (dt[i] == dt_min) { - printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[i], ind[i]); + if (dt[j] == dt_min) { + printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[j], ind[j]); fflush(stdout); } } - } /* i */ + } /* j */ if (config->live_smbh_count > 0) { double min_dt = *std::min_element(dt, dt+N); @@ -1224,10 +1211,10 @@ int main(int argc, char *argv[]) /* load the new values for particles to the local GRAPE's */ /* load the nj particles to the G6 */ - for (int j=0; jlive_smbh_count == 2) { if (config->live_smbh_custom_eps >= 0) { - adjust_bsmbh_softening(eps, config->live_smbh_custom_eps, i_bh1, i_bh2, m_act, x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new); + adjust_bsmbh_softening(eps, config->live_smbh_custom_eps, i_bh1, i_bh2, m[0], m[1], x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new); } if (config->binary_smbh_pn) { // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk - dt_bh_tmp = dt[0]; - - tmp_i = calc_force_pn_BH(m_bh1, x_act_new[i_bh1], v_act_new[i_bh1], s_bh1, - m_bh2, x_act_new[i_bh2], v_act_new[i_bh2], s_bh2, - C_NB, dt_bh_tmp, usedOrNot, + tmp_i = calc_force_pn_BH(m[0], x_act_new[i_bh1], v_act_new[i_bh1], config->smbh1_spin.data(), + m[1], x_act_new[i_bh2], v_act_new[i_bh2], config->smbh2_spin.data(), + config->pn_c, dt[i_bh1], config->pn_usage.data(), (double(*)[3])a_pn1, (double(*)[3])adot_pn1, (double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank); @@ -1447,6 +1433,7 @@ int main(int argc, char *argv[]) get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); #endif + double min_dt = dt_max; for (int i=0; ilive_smbh_count > 0) && (ind_act[i] < config->live_smbh_count)) eta_curr = eta_bh; @@ -1468,7 +1454,6 @@ int main(int argc, char *argv[]) //TODO the below should be moved to a function if (dt_new < dt_min) dt_tmp = dt_min; - if ((dt_new < dt_tmp) && (dt_new > dt_min)) { power = log(dt_new)/log(2.0) - 1; dt_tmp = pow(2.0, (double)power); // TODO why is this casting needed here? @@ -1479,63 +1464,37 @@ int main(int argc, char *argv[]) } if (config->dt_min_warning && (myRank == 0)) { - if (dt_act[i] == dt_min) { - printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]); + if (dt_tmp == dt_min) { + printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt_tmp, ind_act[i]); fflush(stdout); } } + if (dt_tmp < min_dt) min_dt = dt_tmp; - /* BEGIN copy of everything */ - - dt_act[i] = dt_tmp; - t_act[i] = min_t; - - pot_act[i] = pot_act_new[i]; - - x_act[i] = x_act_new[i]; - v_act[i] = v_act_new[i]; - a_act[i] = a_act_new[i]; - adot_act[i] = adot_act_new[i]; - - /* END copy of everything */ - + int j_act = ind_act[i]; + x[j_act] = x_act_new[i]; + v[j_act] = v_act_new[i]; + t[j_act] = min_t; + dt[j_act] = dt_tmp; + pot[j_act] = pot_act_new[i]; + pot_ext[j_act] = pot_act_ext[i]; // ??? + a[j_act] = a_act_new[i]; + adot[j_act] = adot_act_new[i]; } /* i */ /* define the min. dt over all the act. part. and set it also for the BH... */ if (config->live_smbh_count > 0) { - double min_dt = *std::min_element(dt_act, dt_act+n_act); - if (config->live_smbh_count>=1) dt_act[i_bh1] = min_dt; - if (config->live_smbh_count==2) dt_act[i_bh2] = min_dt; + 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)) { //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(i_bh1, i_bh2, ind_act, m_act, x_act_new, v_act_new, pot_act_new, dt_act, n_act, timesteps, time_cur, config->binary_smbh_influence_radius_factor, inf_event); + binary_smbh_influence_sphere_output(i_bh1, i_bh2, ind_act, n_act, timesteps, time_cur, config->binary_smbh_influence_radius_factor, inf_event, &source_particle_list); } - /* Return back the new coordinates + etc... of active part. to the global data... */ - - for (int i=0; i