Skip to content

Generalized 3BBF and bug fixes #75

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ceoconnor1996 wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Generalized 3BBF and bug fixes #75

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9fe4375
generalized/patched 3bbf + various bug fixes
ceoconnor1996 Mar 27, 2026
f1fdeb0
restored mistakenly deleted lines in cmc_dynamics_helper.c
ceoconnor1996 Mar 27, 2026
081e339
Merge branch 'ClusterMonteCarlo:master' into master
ceoconnor1996 Mar 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 15 additions & 5 deletions include/cmc/cmc.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -680,11 +680,16 @@ typedef struct{
#define PARAMDOC_STREAMS "to run the serial version with the given number of random streams - primarily used to mimic the parallel version running with the same no.of processors"
int STREAMS;
/* Meagan - 3bb */
#define PARAMDOC_THREEBODYBINARIES "toggles three-body binary formation (0=off, 1=on)"
#define PARAMDOC_THREEBODYBINARIES "toggles three-body binary formation (0=off, 1=on w/ pairing by mass, 2=on w/ pairing at random)"
/**
* @brief toggles three-body binary formation (0=off, 1=on)
* @brief toggles three-body binary formation (0=off, 1=on w/ pairing by mass, 2=on w/ pairing at random)
*/
int THREEBODYBINARIES;
#define PARAMDOC_BINARY_HARDNESS_POWER "power law index for hardness of newly formed three-body binaries"
/**
* @brief Power-law index for hardness of three-body binaries (default is -3.0).
*/
double BINARY_HARDNESS_POWER;
Comment on lines +688 to +692
#define PARAMDOC_MIN_BINARY_HARDNESS "minimum hardness for newly formed three-body binaries"
/**
* @brief minimum hardness for newly formed three-body binaries
Expand Down Expand Up @@ -1988,7 +1993,7 @@ void binint_do(long k, long kp, double rperi, double w[4], double W, double rcm,

double simul_relax(gsl_rng *rng);
double simul_relax_new(void);
int destroy_bbh(double m1, double m2,double a,double e,double nlocal,double sigma,struct rng_t113_state* rng_st);
int destroy_bbh(double m1,double m2,double a,double e,double nlocal,double sigma,struct rng_t113_state* rng_st);
double simul_relax_new(void);
void calc_sigma_r(long p, long N_LIMIT, double *sig_r, double *sig_sigma, long* sig_n, int r_0_mave_1);
void break_wide_binaries(struct rng_t113_state* rng_st);
Expand All @@ -1998,12 +2003,17 @@ double sigma_r(double r);
// Meagan
/* three-body binary formation */
void sort_three_masses(long sq, long *k1, long *k2, long *k3);
double get_eta(double eta_min, long k1, long k2, long k3, double vrel12[4], double vrel3[4]);
//double get_eta(double eta_min, long k1, long k2, long k3, double vrel12[4], double vrel3[4]);
void calc_3bb_encounter_dyns(long k1, long k2, long k3, double v1[4], double v2[4], double v3[4], double (*vrel12)[4], double (*vrel3)[4], gsl_rng *rng);
void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_binary, double eta_min, double ave_local_mass, double n_local, double sigma_local, double v1[4], double v2[4], double v3[4], double vrel12[4], double vrel3[4], double delta_E_running, gsl_rng *rng);
//void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_binary, double eta_min, double ave_local_mass, double n_local, double sigma_local, double v1[4], double v2[4], double v3[4], double vrel12[4], double vrel3[4], double delta_E_running, gsl_rng *rng);
void calc_sigma_local(long k1, long p, long N_LIMIT, double *ave_local_mass, double *sigma_local);
int remove_old_star(double time, long k);

// Chris O'C -- generalized 3BBF
void make_threebodybinary(double P_3bb, long k1, long k2, long k3, long form_binary, double eta_min, double eta_power, double ave_local_mass, double n_local, double sigma_local, double v1[4], double v2[4], double v3[4], double vrel12[4], double vrel3[4], double delta_E_running, gsl_rng *rng);
void shuffle_three_masses(long sq, long *k1, long *k2, long *k3, gsl_rng *rng);
double get_eta(double eta_min, double eta_power, long k1, long k2, long k3, double vrel12[4], double vrel3[4]);

// Meagan
/* extra output for bhs */
void bh_count(long k);
Expand Down
2 changes: 2 additions & 0 deletions include/cmc/cmc_vars.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,8 @@ _EXTERN_ int BININITKT, STOPATCORECOLLAPSE;
/* Meagan - 3bb */
_EXTERN_ int THREEBODYBINARIES, ONLY_FORM_BH_THREEBODYBINARIES;;
_EXTERN_ double MIN_BINARY_HARDNESS;
// Chris: generalized 3bb
_EXTERN_ double BINARY_HARDNESS_POWER;
_EXTERN_ double BINARY_DISTANCE_BREAKING;
_EXTERN_ int BINARY_BREAKING_MIN;
_EXTERN_ int USE_TT_FILE;
Expand Down
26 changes: 15 additions & 11 deletions src/cmc/cmc_dynamics.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ void dynamics_apply(double dt, gsl_rng *rng)
long sq, k1, k2, k3, form_binary;
double n_threshold, triplet_count, num_triplets_averaged=200;
double ave_local_mass, sigma_local, vrel_ave, v1[4], v2[4], v3[4], vrel12[4], vrel3[4];
double eta_min=MIN_BINARY_HARDNESS, Y1, rate_3bb, rate_ave=0.0, P_3bb, P_ave=0.0;
double eta_min=MIN_BINARY_HARDNESS, eta_power=BINARY_HARDNESS_POWER, Y1, rate_3bb, rate_ave=0.0, P_3bb, P_ave=0.0;
double clight10o7;
double collisions_multiple;

Expand Down Expand Up @@ -66,8 +66,6 @@ void dynamics_apply(double dt, gsl_rng *rng)

pararootfprintf(logfile, "%s(): performing interactions:", __FUNCTION__);



/* Added by Meagan 3/10/11 */
//================================
// Three-body binary formation |
Expand Down Expand Up @@ -96,8 +94,13 @@ void dynamics_apply(double dt, gsl_rng *rng)
{
dt = SaveDt;
form_binary = 0; // reset this to zero; later we decide whether to form a binary, and if so, set form_binary=1
// Sort stars by mass (k1 is most massive)
sort_three_masses(sq, &k1, &k2, &k3);
// Sort stars by mass (k1 is most massive); >Feb 2026 only if using Morscher prescription for 3BBF
if (THREEBODYBINARIES == 1) {
sort_three_masses(sq, &k1, &k2, &k3);
} else {
shuffle_three_masses(sq, &k1, &k2, &k3, rng);
}
//sort_three_masses(sq, &k1, &k2, &k3);
n_local = calc_n_local(get_global_idx(k1), BH_AVEKERNEL, N_LIMIT);
// If density above threshold, check for 3bb formation
if (n_local > n_threshold) {
Expand All @@ -122,17 +125,18 @@ void dynamics_apply(double dt, gsl_rng *rng)
// Calculate RATE of binary formation

// Below is rate_3bb with all the velocity terms vrel_3 and vrel_12 replaced with the averaged local relative velocity, vrel_ave. We did this because we were finding that when we used the actual relative velocities, vrel12, if too large, the 3bb rate would be extremely low and binaries would not form (since it depends strongly on v: v^-9). When we replaced vrel12 with the average relative velocity (over 20 stars), the 3bb formation rate was high enough that binaries would form. We decided to use the average relative velocity for all relative velocity terms, vrel_3 and vrel_12.

// *Note* Factor of 0.5 in front of rate_3bb ensures that our sampling method produces the correct overall analytic 3bb rate
rate_3bb = 0.5 * sqrt(2) * sqr(PI) * sqr(n_local) * pow(vrel_ave, -9) * pow(((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)]) * madhoc), 5.0) * pow(eta_min, -5.5) * (1.0 + 2.0*eta_min) * (1.0 + 2.0 * ((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)] + star_m[get_global_idx(k3)]) / (star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)])) * eta_min);
// old rate from Morscher et al. 2013
//rate_3bb = 0.5 * sqrt(2) * sqr(PI) * sqr(n_local) * pow(vrel_ave, -9) * pow(((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)]) * madhoc), 5.0) * pow(eta_min, -eta_power) * (1.0 + 2.0*eta_min) * (1.0 + 2.0 * ((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)] + star_m[get_global_idx(k3)]) / (star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)])) * eta_min);
// updated rate from Atallah et al. 2024, plus user-specified power law
rate_3bb = 0.1328 * sqrt(10.) * sqr(n_local) * pow(vrel_ave,-9) * pow(((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)]) * madhoc), 5.0) * pow(eta_min, eta_power) * (1.0 + 2.0*eta_min) * (1.0 + 2.0 * ((star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)] + star_m[get_global_idx(k3)]) / (star_m[get_global_idx(k1)] + star_m[get_global_idx(k2)])) * eta_min);

// Calculate PROBABILITY of binary formation
P_3bb = rate_3bb * (dt * ((double) clus.N_STAR)/log(GAMMA*((double) clus.N_STAR)));

// print info on probability calculation, for innermost 200 triplets, in each timestep
if (triplet_count <= num_triplets_averaged) {
P_ave += P_3bb;
rate_ave += rate_3bb;
P_ave += P_3bb;
rate_ave += rate_3bb;
}
/* //MPI: Commenting out since this will cause a deadlock, and is also this needs a reduction for each star.
// For each timestep, output average rate and probability of 3bb formation,
Expand Down Expand Up @@ -171,7 +175,7 @@ void dynamics_apply(double dt, gsl_rng *rng)
/* Here is where new binary properties are calculated.
if form_binary=1, the new binary will be created, but if form_binary=0,
the interaction will be logged in lightcollision log */
make_threebodybinary(P_3bb, k1, k2, k3, form_binary, eta_min, ave_local_mass, n_local, sigma_local, v1, v2, v3, vrel12, vrel3, delta_E_3bb, rng);
make_threebodybinary(P_3bb, k1, k2, k3, form_binary, eta_min, eta_power, ave_local_mass, n_local, sigma_local, v1, v2, v3, vrel12, vrel3, delta_E_3bb, rng);
} else { // Probability of 3bb formation too low ==> No binary formed
form_binary = 0;
/* do nothing; check 3bb formation for next 3 stars */
Expand Down
Loading
Loading