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

Moved SMBH de-softening to a function (still not pretty) and improved init script

Browse source
This commit is contained in:
Yohai Meiron 2020-04-05 19:54:10 -04:00
parent 56abe820c3
commit d9e3aea243
2 changed files with 100 additions and 47 deletions
Download patch file Download diff file Expand all files Collapse all files

49
init.py
View file

@ -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)

98
phigrape.cpp
View file

@ -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; i<n_act; i++) {
int j_act = ind_act[i];
// NOTICE much of these are just not needed!
m_act[i] = m[j_act];
x_act[i] = x[j_act];
v_act[i] = v[j_act];
@ -1370,48 +1419,7 @@ int main(int argc, char *argv[])
if (config->live_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... */