Enabled triple (plus) SMBHs
This commit is contained in:
parent
fcdafd94a6
commit
633c82b917
4 changed files with 75 additions and 148 deletions
141
black_holes.cpp
141
black_holes.cpp
|
|
@ -10,64 +10,39 @@ double L[3]; // needed in pn_bh_spin.c
|
|||
#undef SQR
|
||||
/* END legacy inclusion */
|
||||
|
||||
void two_body_gravity(
|
||||
const double m1, const double3 x1, const double3 v1,
|
||||
const double m2, const double3 x2, const double3 v2,
|
||||
const double eps,
|
||||
double& pot1, double3& acc1, double3& jrk1,
|
||||
double& pot2, double3& acc2, double3& jrk2)
|
||||
void two_body_gravity(const Particle_ref& j, const Particle_ref& i, const double eps, double& pot, double3& acc, double3& jrk)
|
||||
{
|
||||
double3 dx = x1 - x2;
|
||||
double3 dv = v1 - v2;
|
||||
double3 dx = i.x - j.x;
|
||||
double3 dv = i.v - j.v;
|
||||
double r2 = dx.norm2() + eps*eps;
|
||||
double r = sqrt(r2);
|
||||
double r3 = r2*r;
|
||||
double r4 = r2*r2;
|
||||
double RP = 3*(dx*dv)/r;
|
||||
|
||||
pot1 = -m2/r;
|
||||
pot2 = -m1/r;
|
||||
acc1 = -m2*dx/r3;
|
||||
acc2 = m1*dx/r3;
|
||||
jrk1 = -m2*(dv/r3 - RP*dx/r4);
|
||||
jrk2 = m1*(dv/r3 - RP*dx/r4);
|
||||
pot = -j.m/r;
|
||||
acc = -j.m*dx/r3;
|
||||
jrk = -j.m*(dv/r3 - RP*dx/r4);
|
||||
}
|
||||
|
||||
void Black_hole_physics::adjust_softening(
|
||||
double& pot1, double& pot2,
|
||||
double3& acc1, double3& acc2,
|
||||
double3& jrk1, double3& jrk2)
|
||||
void Black_hole_physics::adjust_softening(const std::vector<Particle_ref>& particles)
|
||||
{
|
||||
if (eps_new < 0) return;
|
||||
// calculate and "minus" the BH <-> BH softened pot, acc & jerk
|
||||
two_body_gravity(
|
||||
m1, x1, v1,
|
||||
m2, x2, v2,
|
||||
eps_old,
|
||||
bbh_grav.pot1, bbh_grav.a1, bbh_grav.adot1,
|
||||
bbh_grav.pot2, bbh_grav.a2, bbh_grav.adot2);
|
||||
pot1 -= bbh_grav.pot1;
|
||||
pot2 -= bbh_grav.pot2;
|
||||
acc1 -= bbh_grav.a1;
|
||||
acc2 -= bbh_grav.a2;
|
||||
jrk1 -= bbh_grav.adot1;
|
||||
jrk2 -= bbh_grav.adot2;
|
||||
|
||||
// calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk
|
||||
two_body_gravity(
|
||||
m1, x1, v1,
|
||||
m2, x2, v2,
|
||||
eps_new,
|
||||
bbh_grav.pot1, bbh_grav.a1, bbh_grav.adot1,
|
||||
bbh_grav.pot2, bbh_grav.a2, bbh_grav.adot2);
|
||||
pot1 += bbh_grav.pot1;
|
||||
pot2 += bbh_grav.pot2;
|
||||
acc1 += bbh_grav.a1;
|
||||
acc2 += bbh_grav.a2;
|
||||
jrk1 += bbh_grav.adot1;
|
||||
jrk2 += bbh_grav.adot2;
|
||||
double pot_old, pot_new;
|
||||
double3 acc_old, acc_new, jrk_old, jrk_new;
|
||||
for (int j = 0; j < particles.size(); j++) {
|
||||
for (int i = 0; i < particles.size(); i++) {
|
||||
if (i == j) continue;
|
||||
two_body_gravity(particles[j], particles[i], eps_old, pot_old, acc_old, jrk_old);
|
||||
two_body_gravity(particles[j], particles[i], eps_new, pot_new, acc_new, jrk_new);
|
||||
particles[i].pot += pot_new - pot_old;
|
||||
particles[i].acc += acc_new - acc_old;
|
||||
particles[i].jrk += jrk_new - jrk_old;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#if 0
|
||||
void Black_hole_physics::adjust_post_newtonian(
|
||||
const double dt_bh, // pn_usage should be const
|
||||
double3& acc1, double3& acc2,
|
||||
|
|
@ -76,8 +51,8 @@ void Black_hole_physics::adjust_post_newtonian(
|
|||
// calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
|
||||
// TODO maybe have the PN terms as local variables here?
|
||||
int tmp;
|
||||
tmp = calc_force_pn_BH(m1, x1, v1, bbh_grav.spin1,
|
||||
m2, x2, v2, bbh_grav.spin2,
|
||||
tmp = calc_force_pn_BH(masses[0], x1, v1, bbh_grav.spin1,
|
||||
masses[1], x2, v2, bbh_grav.spin2,
|
||||
c, dt_bh, pn_usage,
|
||||
(double(*)[3])bbh_grav.a_pn1, (double(*)[3])bbh_grav.adot_pn1,
|
||||
(double(*)[3])bbh_grav.a_pn2, (double(*)[3])bbh_grav.adot_pn2, myRank, rootRank);
|
||||
|
|
@ -98,72 +73,22 @@ void Black_hole_physics::adjust_post_newtonian(
|
|||
jrk1 += jrk1_corr;
|
||||
jrk2 += jrk2_corr;
|
||||
}
|
||||
#endif
|
||||
|
||||
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)
|
||||
void Black_hole_physics::write_bh_data(const double time_cur, const int count, 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
|
||||
// (the arguments to this function) and optionall internal data to this
|
||||
// object (the most recently calculated force and post-Newtonian terms).
|
||||
auto out = fopen("bh.dat", "a");
|
||||
if (count == 2) {
|
||||
for (int i=0; i < 2; i++) {
|
||||
double3 *a_pn, *adot_pn, a_bh, adot_bh;
|
||||
double pot_bh;
|
||||
if (i==0) {
|
||||
a_pn = bbh_grav.a_pn1;
|
||||
adot_pn = bbh_grav.adot_pn1;
|
||||
pot_bh = bbh_grav.pot1;
|
||||
a_bh = bbh_grav.a1;
|
||||
adot_bh = bbh_grav.adot1;
|
||||
} else {
|
||||
a_pn = bbh_grav.a_pn2;
|
||||
adot_pn = bbh_grav.adot_pn2;
|
||||
pot_bh = bbh_grav.pot2;
|
||||
a_bh = bbh_grav.a2;
|
||||
adot_bh = bbh_grav.adot2;
|
||||
}
|
||||
if (eps_new >= 0) {
|
||||
fprintf(out, "%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E ",
|
||||
time_cur, m[i],
|
||||
x[i][0], x[i][1], x[i][2], x[i].norm(),
|
||||
v[i][0], v[i][1], v[i][2], v[i].norm(),
|
||||
pot[i],
|
||||
a[i][0], a[i][1], a[i][2], a[i].norm(),
|
||||
adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
|
||||
dt[i],
|
||||
pot_bh,
|
||||
a_bh[0], a_bh[1], a_bh[2], a_bh.norm(),
|
||||
adot_bh[0], adot_bh[1], adot_bh[2], adot_bh.norm());
|
||||
if (c > 0) {
|
||||
fprintf(out, "\t");
|
||||
for (int pn_idx=0; pn_idx < 7; pn_idx++) {
|
||||
fprintf(out, "\t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E ", a_pn[pn_idx][0], a_pn[pn_idx][1], a_pn[pn_idx][2], a_pn[pn_idx].norm(), adot_pn[pn_idx][0], adot_pn[pn_idx][1], adot_pn[pn_idx][2], adot_pn[pn_idx].norm());
|
||||
}
|
||||
}
|
||||
fprintf(out, "\n");
|
||||
} else {
|
||||
fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \n",
|
||||
time_cur, m[i],
|
||||
x[i][0], x[i][1], x[i][2], x[i].norm(),
|
||||
v[i][0], v[i][1], v[i][2], v[i].norm(),
|
||||
pot[i],
|
||||
a[i][0], a[i][1], a[i][2], a[i].norm(),
|
||||
adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
|
||||
dt[i]);
|
||||
}
|
||||
}
|
||||
fprintf(out, "\n");
|
||||
} else if (count == 1) {
|
||||
for (int i = 0; i < count; i++) {
|
||||
fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \n",
|
||||
time_cur, m[0],
|
||||
x[0][0], x[0][1], x[0][2], x[0].norm(),
|
||||
v[0][0], v[0][1], v[0][2], v[0].norm(),
|
||||
pot[0],
|
||||
a[0][0], a[0][1], a[0][2], a[0].norm(),
|
||||
adot[0][0], adot[0][1], adot[0][2], adot[0].norm(),
|
||||
dt[0]);
|
||||
fprintf(out,"\n");
|
||||
time_cur, m[i],
|
||||
x[i][0], x[i][1], x[i][2], x[i].norm(),
|
||||
v[i][0], v[i][1], v[i][2], v[i].norm(),
|
||||
pot[i],
|
||||
a[i][0], a[i][1], a[i][2], a[i].norm(),
|
||||
adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
|
||||
dt[i]);
|
||||
}
|
||||
fprintf(out, "\n");
|
||||
fclose(out);
|
||||
}
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue