From d9e3aea24397b827f53878a381eaa94184ab0fb4 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Sun, 5 Apr 2020 19:54:10 -0400 Subject: [PATCH] Moved SMBH de-softening to a function (still not pretty) and improved init script --- init.py | 49 +++++++++++++++++++++++++- phigrape.cpp | 98 ++++++++++++++++++++++++++++------------------------ 2 files changed, 100 insertions(+), 47 deletions(-) diff --git a/init.py b/init.py index e5aa76c..a12d55c 100644 --- a/init.py +++ b/init.py @@ -33,6 +33,38 @@ def gen_plum(N, seed=None, RMAX=10): i += 1 return particle_list +def kepler_to_cartesian(a, e, i, Omega, w, nu, G=1.0, M=1.0): + def to_arrays(*args): + 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) + 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)] + Q = [-np.sin(w)*np.cos(Omega) - np.cos(w)*np.cos(i)*np.sin(Omega), + -np.sin(w)*np.sin(Omega) + np.cos(w)*np.cos(i)*np.cos(Omega), + np.sin(i)*np.cos(w)] + cosnu = np.cos(nu) + cosE = (e+cosnu)/(1+e*cosnu) + E = np.arccos(cosE) + E[nu > np.pi] = 2*np.pi - E[nu > np.pi] + X = a*((np.cos(E)-e)*P + np.sqrt(1-e**2)*np.sin(E)*Q) + V = (np.sqrt(G*M/a)/(1-e*np.cos(E)))*(-np.sin(E)*P + np.sqrt(1-e**2)*np.cos(E)*Q) + if X.shape[1]==1: + X=X[:,0] + V=V[:,0] + return X.T, V.T + +def generate_binary(a, e, i, Omega, w, nu, m1, m2): + X, V = kepler_to_cartesian(a, e, i, Omega, w, nu, M=m1+m2) + q = np.double(m2)/np.double(m1) + X1 = -q/(q+1)*X + V1 = -q/(q+1)*V + X2 = 1/(q+1)*X + V2 = 1/(q+1)*V + return X1, V1, X2, V2 + def write_phi_grape_config(**kargs): if 'file_name' in kargs: file_name = kargs['file_name'] else: file_name = 'phi-GRAPE.cfg' @@ -86,6 +118,7 @@ if __name__=='__main__': parser.add_argument('--dt_bh', type=np.double, default=.125, help='interval for BH output (bh.dat & bh_nb.dat)') parser.add_argument('--eta', type=np.double, default=.01, help='parameter for timestep determination (0.02 or 0.01)') parser.add_argument('--frac', type=np.double, default=0, help='fraction of collisional particles (by angular momentum)') + parser.add_argument('--bsmbh', type=bool, default=0, help='generate a binary supermassive black hole (parameters hardcoded in the script)') args = parser.parse_args() try: @@ -97,9 +130,23 @@ if __name__=='__main__': write_phi_grape_config(**vars(args)) particle_list = gen_plum(N, seed=args.seed) + m = np.ones(N)/N - write_particles(particle_list) + if args.bsmbh: + m1, m2 = 0.075, 0.025 + a, e, i, Omega, w, nu = 0.001, 0.5, 0, 0, 0, 0 + X1, V1, X2, V2 = generate_binary(a, e, i, Omega, w, nu, m1, m2) + m[:2] = m1, m2 + m[2:] = 1/(N-2) + particle_list[0,:3] = X1 + particle_list[0,3:] = V1 + particle_list[1,:3] = X2 + particle_list[1,3:] = V2 + + write_particles(particle_list, m=m) mask = gen_mask(particle_list, args.frac) + if args.bsmbh: + mask[:2] = 3 write_mask(mask) diff --git a/phigrape.cpp b/phigrape.cpp index 090bcfa..dfe3e13 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -666,6 +666,54 @@ void binary_smbh_influence_sphere_output(int i_bh1, int i_bh2, int ind_act[], do } +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[]) +{ + 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]; + + double3 x_bh2 = x_act_new[i_bh2]; + double3 v_bh2 = v_act_new[i_bh2]; + + // calculate and "minus" the BH <-> BH softened pot, acc & jerk + + double tmp_i; + tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1, + m_bh2, x_bh2, v_bh2, + eps, + &pot_bh1, a_bh1, adot_bh1, + &pot_bh2, a_bh2, adot_bh2); + + pot_act_new[i_bh1] -= pot_bh1; + pot_act_new[i_bh2] -= pot_bh2; + + a_act_new[i_bh1] -= a_bh1; + a_act_new[i_bh2] -= a_bh2; + + adot_act_new[i_bh1] -= adot_bh1; + adot_act_new[i_bh2] -= adot_bh2; + + // calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk + + tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1, + m_bh2, x_bh2, v_bh2, + eps_bh, + &pot_bh1, a_bh1, adot_bh1, + &pot_bh2, a_bh2, adot_bh2); + + pot_act_new[i_bh1] += pot_bh1; + pot_act_new[i_bh2] += pot_bh2; + + a_act_new[i_bh1] += a_bh1; + a_act_new[i_bh2] += a_bh2; + + adot_act_new[i_bh1] += adot_bh1; + adot_act_new[i_bh2] += adot_bh2; + +} + int main(int argc, char *argv[]) { int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank, @@ -1306,6 +1354,7 @@ int main(int argc, char *argv[]) #endif for (int i=0; ilive_smbh_count == 2) { if (config->live_smbh_custom_eps >= 0) { - m_bh1 = m_act[i_bh1]; - m_bh2 = m_act[i_bh2]; - - x_bh1 = x_act_new[i_bh1]; - v_bh1 = v_act_new[i_bh1]; - - x_bh2 = x_act_new[i_bh2]; - v_bh2 = v_act_new[i_bh2]; - - // calculate and "minus" the BH <-> BH softened pot, acc & jerk - - tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1, - m_bh2, x_bh2, v_bh2, - eps, - &pot_bh1, a_bh1, adot_bh1, - &pot_bh2, a_bh2, adot_bh2); - - pot_act_new[i_bh1] -= pot_bh1; - pot_act_new[i_bh2] -= pot_bh2; - - a_act_new[i_bh1] -= a_bh1; - a_act_new[i_bh2] -= a_bh2; - - adot_act_new[i_bh1] -= adot_bh1; - adot_act_new[i_bh2] -= adot_bh2; - - // calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk - - tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1, - m_bh2, x_bh2, v_bh2, - config->live_smbh_custom_eps, - &pot_bh1, a_bh1, adot_bh1, - &pot_bh2, a_bh2, adot_bh2); - - pot_act_new[i_bh1] += pot_bh1; - pot_act_new[i_bh2] += pot_bh2; - - a_act_new[i_bh1] += a_bh1; - a_act_new[i_bh2] += a_bh2; - - adot_act_new[i_bh1] += adot_bh1; - adot_act_new[i_bh2] += adot_bh2; + 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); } if (config->binary_smbh_pn) { @@ -1498,11 +1506,9 @@ int main(int argc, char *argv[]) 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 */ - - } /* i */ /* define the min. dt over all the act. part. and set it also for the BH... */ @@ -1515,7 +1521,7 @@ int main(int argc, char *argv[]) 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, v_act, pot_act, 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, 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); } /* Return back the new coordinates + etc... of active part. to the global data... */