Commit 32e7c652 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Also apply the distance criterion when checking for BH mergers. Read the...

Also apply the distance criterion when checking for BH mergers. Read the correct value of G from the physical constant structure.
parent 6ce26301
......@@ -358,12 +358,12 @@ int main(int argc, char *argv[]) {
return 1;
}
if (with_black_holes && !with_external_gravity && !with_self_gravity) {
if (with_black_holes && !with_self_gravity) {
if (myrank == 0) {
argparse_usage(&argparse);
printf(
"\nError: Cannot process black holes without gravity, -g or -G "
"must be chosen.\n");
"\nError: Cannot process black holes without self-gravity, -G must "
"be chosen.\n");
}
return 1;
}
......@@ -799,11 +799,10 @@ int main(int argc, char *argv[]) {
bzero(&black_holes_properties, sizeof(struct black_holes_props));
/* Initialise the gravity properties */
bzero(&gravity_properties, sizeof(struct gravity_props));
if (with_self_gravity)
gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology,
periodic);
else
bzero(&gravity_properties, sizeof(struct gravity_props));
gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
with_cosmology, periodic);
/* Initialise the cooling function properties */
#ifdef COOLING_NONE
......
......@@ -410,8 +410,9 @@ int main(int argc, char *argv[]) {
cosmology_init(params, &us, &prog_const, &cosmo);
/* Initialise the gravity scheme */
gravity_props_init(&gravity_properties, params, &cosmo, /*with_cosmology=*/1,
periodic);
bzero(&gravity_properties, sizeof(struct gravity_props));
gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
/*with_cosmology=*/1, periodic);
/* Initialise the FOF properties */
bzero(&fof_properties, sizeof(struct fof_props));
......
......@@ -112,6 +112,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const struct bpart *restrict bi,
struct bpart *restrict bj,
const struct cosmology *cosmo,
const struct gravity_props *grav_props,
const integertime_t ti_current) {}
/**
......
......@@ -196,6 +196,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
const struct bpart *restrict bi,
struct bpart *restrict bj,
const struct cosmology *cosmo,
const struct gravity_props *grav_props,
const integertime_t ti_current) {
/* Compute relative peculiar velocity between the two BHs
......@@ -215,7 +216,11 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
h = hj;
}
const float G = 43.; // MATTHIEU: Fix this!!!
const float max_distance2 = kernel_gravity_softening_plummer_equivalent_inv *
kernel_gravity_softening_plummer_equivalent_inv *
9.f * grav_props->epsilon_cur2;
const float G_Newton = grav_props->G_Newton;
/* The BH with the smaller mass will be merged onto the one with the
* larger mass.
......@@ -224,12 +229,10 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
if ((bj->subgrid_mass < bi->subgrid_mass) ||
(bj->subgrid_mass == bi->subgrid_mass && bj->id < bi->id)) {
/* Merge if gravitationally bound
/* Merge if gravitationally bound AND if within max distance
* Note that we use the kernel support here as the size and not just the
* smoothing length */
if (v2_pec < G * M / (kernel_gamma * h)) {
// MATTHIEU: Also add distance check: factor * Plummer softening
if (v2_pec < G_Newton * M / (kernel_gamma * h) && (r2 < max_distance2)) {
/* This particle is swallowed by the BH with the largest ID of all the
* candidates wanting to swallow it */
......
......@@ -39,6 +39,7 @@
#define gravity_props_default_rebuild_frequency 0.01f
void gravity_props_init(struct gravity_props *p, struct swift_params *params,
const struct phys_const *phys_const,
const struct cosmology *cosmo, int with_cosmology,
int periodic) {
......@@ -98,6 +99,9 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
p->epsilon_comoving = p->epsilon_max_physical;
}
/* Copy over the gravitational constant */
p->G_Newton = phys_const->const_newton_G;
/* Set the softening to the current time */
gravity_props_update(p, cosmo);
}
......
......@@ -31,6 +31,7 @@
/* Forward declarations */
struct cosmology;
struct phys_const;
struct swift_params;
/**
......@@ -84,10 +85,15 @@ struct gravity_props {
/*! Cube of the inverse of current oftening length */
float epsilon_cur_inv3;
/*! Gravitational constant (in internal units, copied from the physical
* constants) */
float G_Newton;
};
void gravity_props_print(const struct gravity_props *p);
void gravity_props_init(struct gravity_props *p, struct swift_params *params,
const struct phys_const *phys_const,
const struct cosmology *cosmo, int with_cosmology,
int periodic);
void gravity_props_update(struct gravity_props *p,
......
......@@ -209,7 +209,8 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
#endif
if (r2 < hig2) {
IACT_BH_BH(r2, dx, hi, hj, bi, bj, cosmo, ti_current);
IACT_BH_BH(r2, dx, hi, hj, bi, bj, cosmo, e->gravity_properties,
ti_current);
}
} /* loop over the bparts in ci. */
} /* loop over the bparts in ci. */
......@@ -355,7 +356,8 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
#endif
if (r2 < hig2) {
IACT_BH_BH(r2, dx, hi, hj, bi, bj, cosmo, ti_current);
IACT_BH_BH(r2, dx, hi, hj, bi, bj, cosmo, e->gravity_properties,
ti_current);
}
} /* loop over the bparts in cj. */
} /* loop over the bparts in ci. */
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment