From e9455c0a0e62396cf78103547e3a3b3945837946 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Wed, 4 Mar 2020 13:46:04 -0500 Subject: [PATCH] Started cleanup process --- Makefile | 14 +- phi-GRAPE.c | 7293 -------------------------------------------------- phigrape.cpp | 6934 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 6941 insertions(+), 7300 deletions(-) delete mode 100644 phi-GRAPE.c create mode 100644 phigrape.cpp diff --git a/Makefile b/Makefile index f8fe666..20340bf 100644 --- a/Makefile +++ b/Makefile @@ -13,18 +13,18 @@ yebisu: GRAPEHOME = ../yebisu yebisu: GRAPELIB = -L$(GRAPEHOME) -lyebisug6 GRAPEINC = -I$(GRAPEHOME) -CFLAGS ?= -mcmodel=large -CFLAGS += -O$(OPTIMIZATION) +CXXFLAGS ?= -mcmodel=medium +CFFFLAGS += -O$(OPTIMIZATION) INC = $(GRAPEINC) $(CUDAINC) -LIB = $(GRAPELIB) $(CUDALIB) -lm -lgcc -lgfortran -lstdc++ -MPICC ?= mpicc -EXECUTABLE ?= phi-GRAPE.exe +LIB = $(GRAPELIB) $(CUDALIB) -lm +MPICXX ?= mpic++ +EXECUTABLE ?= phigrape default: - $(MPICC) $(CPPFLAGS) $(CFLAGS) -DETICS_DTSCF=$(ETICS_DTSCF) $(INC) phi-GRAPE.c -o $(EXECUTABLE) $(LIB) + $(MPICXX) $(CPPFLAGS) $(CXXFLAGS) -DETICS_DTSCF=$(ETICS_DTSCF) $(INC) phigrape.cpp -o $(EXECUTABLE) $(LIB) yebisu: CPPFLAGS := $(filter-out -DETICS, $(CPPFLAGS)) yebisu: default clean: - rm -f *.o phi-GRAPE.exe + rm -f *.o phigrape diff --git a/phi-GRAPE.c b/phi-GRAPE.c deleted file mode 100644 index 236ee7e..0000000 --- a/phi-GRAPE.c +++ /dev/null @@ -1,7293 +0,0 @@ -/***************************************************************************** - File Name : "phi-GRAPE/GPU.c" // BH (1 || 2) + ACC + EJECT - : - Contents : N-body code with integration by individual block time step - : together with the parallel using of GRAPE6a board's. - : - : Added the GPU support via SAPPORO library. - : - : Normalization to the physical units!!! - : - : External Potential added - : Plummer-Kuzmin: Bulge, Disk, Halo - : Kharchenko+Andreas... - : - : SC extra POT for Bek SC test runs... - : - : Rebuced to the Single BH -> Plummer - : Andreas+Fazeel... - : - : Stellar evolution added - : Stellar lifetimes: Raiteri, Villata & Navarro (1996) - : IMS mass loss: van den Hoeg & Groenewegen (1997) - : - : STARDESTR_EXT: Tidal disruption of stars by external BH... - : Chingis, Denis & Maxim... - : - : STARDESTR: Tidal disruption of stars by BH... - : Jose, Li Shuo & Shiyan Zhong - : - : STARDISK: Drag force... - : Chingis, Denis & Maxim... - : - : STARDISK: variable hz = HZ*(R/R_crit) up to R_crit... - : Taras, Andreas... - : - : Live BH (1 || 2) + ACC + EJECT... - : Li Shuo & Shiyan Zhong - : - : dt_min for BH (1 || 2)... - : - : added the PN calculus for the BBH - : PN0, PN1, PN2, PN2.5 (coded on the base of - : Gabor Kupi original routine) - : - : added the "name" array... - : - : added the GMC's calculus (GMC on CPU; GMC2 on GPU) - : for Alexey SC runs... and also for Fazeel Zurich runs... - : - : CPU_TIMELIMIT added for the Julich MW cluster runs... - : - Coded by : Peter Berczik - Version number : 19.04 - Last redaction : 2019.04.16 12:55 -*****************************************************************************/ -//#define CPU_RUNLIMIT 28000 // 24h = 24*60*60 = 86400 -#define CPU_RUNLIMIT 255000 // 3 days = 3*24*60*60 = 259200 -//#define CPU_RUNLIMIT 100 // 9.5 min -//#define CPU_RUNLIMIT 100 // 9.5 min - -//#define GMC // add the GMC's to the run -//#define GMC2 // add the GMC's to the run - -//#define NORM // Physical normalization -//#define STEVOL // Stellar evolution Do not tuch! Defined inside the Makefiles !!! -//#define STEVOL_SSE // Stellar evolution with SSE Do not tuch! Defined inside the Makefiles !!! - -//#define ADD_BH1 // add the Single BH - -#define ADD_BH2 // add the Binary BH's -#define ADD_N_BH // eps_BH = 0.0, but added only the Newtonian forces -// #define ADD_PN_BH // extra - added also the Post-Newton forces -// #define ADD_SPIN_BH // extra - added the SPIN for the BH's - DEFAULT !!! - -// #define BH_OUT // extra output for BH's (live) -// #define BH_OUT_NB // extra output for the BH's neighbours (live) - -// #define BBH_INF // BBH influence sphere... -// #define R_INF 10.0 // Factor for the influence sphere... if( R < R_INF * DR_BBH ) -// #define R_INF2 (R_INF*R_INF) - -// #define DT_MIN_WARNING // dt < dt_min warning !!! - -//#define BH_OUT_NB_EXT // extra output for the BH's neighbours (external) - -//#define SAPPORO // Gaburov. Do not tuch! Defined inside the Makefiles !!! -//#define YEBISU // Nitadori. Do not tuch! Defined inside the Makefiles !!! -//#define GUINNESS // Nakasato. Do not tuch! Defined inside the Makefiles !!! - -//#define STARDESTR // disruption of stars by BH tidal forces -//#define R_TIDAL 1.0E-03 // Tidal radius of the BH's -//#define RMAX_OUT // extra data for def. r_max... - -//#define STARDESTR_EXT // disruption of stars by external BH tidal forces -//#define R_0 0.22 // Outer cut of the disk -//#define R_ACCR 0.0050 // Tidal Accr. rad of the BH's in units of R_0 -//#define R_TIDAL (R_0*R_ACCR) // Tidal radius of the BH's in NB units - -//#define STARDISK // analytic gas disk around BH + drag -//#define HZ (0.001*R_0) // Disk thickness... in NB units -//#define R_CRIT 0.0257314 // Critical radius of the disk (vert SG = vert BH force) -//#define R_CRIT 0.0 // i.e. hz=HZ=const... - -//#define SPH_GAS // SPH gas disk around BH + drag; experimental stuff !!! - -//#define EXTPOT // external potential (BH? or galactic?) -//#define EXTPOT_BH // BH - usually NB units - -//#define EXTPOT_GAL // Galactic B+D+H PK - usually physical units - -//#define EXTPOT_GAL_DEH // Dehnen Galactic - usually physical units -//#define EXTPOT_GAL_LOG // Log Galactic - usually physical units - -//#define EXTPOT_SC // SC extra POT for Bek test runs... -//#define DATAFILE eff0.05.tsf10.00.tab -//#define M_SC_DIM 100001 - -//#define CMCORR // CM correction in the zero step and in every dt_contr - -//#define DEBUG // debug information to files & screen -// #define DEBUG_extra // extra debug information to files & screen -//#define LAGR_RAD // write out lagr-rad.dat - -//#define LIMITS // for "0" mass particles !!! -//#define R_LIMITS 1.0E+03 // for "0" mass particles !!! -//#define COORDS 1010.0 // for "0" mass particles !!! - -//#define LIMITS_NEW // for "0" mass particles !!! - -//#define LIMITS_ALL_BEG // for ALL particles at the beginning... -//#define LIMITS_ALL // for ALL particles -//#define R_LIMITS_ALL 1.0E+03 // for ALL particles - -#ifdef ETICS -#include "grapite.h" -// why do we need CEP as a compilaion flag... just have it always on when ETICS is on. IF there is no CEP, there should be a graceful skipping of those operations. -//#define ETICS_CEP -#ifndef ETICS_DTSCF -#error "ETICS_DTSCF must be defined" -#endif -#endif - -#define TIMING - -#define ETA_S_CORR 4.0 -#define ETA_BH_CORR 4.0 - -#define DTMAXPOWER -3.0 -#define DTMINPOWER -36.0 - -/* - -3.0 0.125 - -4.0 0.0625 - -5.0 0.03125 - -7.0 ~1e-2 --10.0 ~1e-3 -............. --20.0 ~1e-6 --23.0 ~1e-7 --26.0 ~1e-8 --30.0 ~1e-9 -*/ - -//#define AMD - -// #define ACT_DEF_LL - -#if defined(ACT_DEF_LL) && defined(ACT_DEF_GRAPITE) -#error "Contradicting preprocessor flags!" -#endif - -/****************************************************************************/ -#include -#include -#include -#include -#include -#include -#include - -/* -double aaa; -double aaapars[5]; - -extern void qwerty_(double *aaa, double *aaapars); -*/ - -#ifdef SAPPORO -#include -#include -#define GPU_PIPE 256 -#else -# ifdef YEBISU -# define G6_NPIPE 2048 -# else -# define G6_NPIPE 48 -# endif -#include "grape6.h" -#endif - -#ifdef GUINNESS -#define GPU_PIPE 128 -#endif - -#include - -/* Some "good" functions and constants... */ -#define SIG(x) ( ((x)<0) ? (-1):(1) ) -#define ABS(x) ( ((x)<0) ? (-x):(x) ) -#define MAX(a,b) ( ((a)>(b)) ? (a):(b) ) -#define MIN(a,b) ( ((a)<(b)) ? (a):(b) ) -#define SQR(x) ( (x)*(x) ) -#define POW3(x) ( (x)*SQR(x) ) - -#define Pi 3.14159265358979323846 -#define TWOPi 6.283185307179 -#define sqrt_TWOPi 2.506628274631 - -#ifdef NORM -//http://pdg.lbl.gov/2015/reviews/rpp2015-rev-astrophysical-constants.pdf - -#define G 6.67388E-11 // (m/s^2) * (m^2/kg) -#define Msol 1.988489E+30 // kg -#define Rsol 6.957E+08 // m -#define AU 149597870700.0 // m -#define pc 3.08567758149E+16 // m -#define Year 31556925.2 // s -#define c_feny 299792458.0 // m/s - -#define kpc (1.0E+03*pc) // m -#define km 1.0E+03 // km -> m -#define cm3 1.0E-06 // cm^3 -> m^3 -#define Myr (1.0E+06*Year) // s -#define Gyr (1.0E+09*Year) // s -#define R_gas 8.31447215 // J/(K*mol) -#define k_gas 1.380650424E-23 // J/K -#define N_A 6.022141510E+23 // 1/mol -#define mu 1.6605388628E-27 // kg -#define mp 1.67262163783E-27 // kg -#define me 9.1093821545E-31 // kg - -#define pc2 (pc*pc) -#define pc3 (pc*pc*pc) -#define kpc2 (kpc*kpc) -#define kpc3 (kpc*kpc*kpc) -#endif - -/* - 1KB = 1024 - 2KB = 2048 - 4KB = 4096 - 8KB = 8192 - 16KB = 16384 - 32KB = 32768 - 64KB = 65536 - 128KB = 131072 - 256KB = 262144 - 512KB = 524288 - 1024KB = 1048576 -> 1MB -*/ - -#define KB 1024 -#define MB (KB*KB) - -#define N_MAX (6*MB) -#define N_MAX_loc (2*MB) - -//#define N_MAX (1200000) -//#define N_MAX_loc (1200000) - -//#define P_MAX 32 - -//#ifdef DEBUG - -/* extra code & data for DEBUG... */ - -#include "debug.h" - -//#ifdef DEBUG_extra -int ind_sort[N_MAX]; -double var_sort[N_MAX]; -//#endif - -//#endif - -#ifdef LAGR_RAD -int lagr_rad_N = 22; -double mass_frac[] = { 0.0001, 0.0003, 0.0005, 0.001, 0.003, 0.005, 0.01, 0.03, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 0.9999, 1.0 }; -double lagr_rad[22]; -double m_tot, m_cum, n_cum; -#endif - -int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank, - i, j, k, ni, nj, diskstep=0, power, jjj, iii, - skip_con=0, tmp_i; - -double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, - dt_bh, t_bh=0.0, dt_bh_tmp, - t_end, time_cur, dt_min, dt_max, min_t, min_t_loc, dt_new, min_dt, - eta_s, eta, eta_bh, - E_pot, E_pot_ext, E_kin, E_tot, E_tot_0, DE_tot, - E_tot_corr, E_tot_corr_0, DE_tot_corr, - E_tot_corr_sd, E_tot_corr_sd_0, DE_tot_corr_sd, - E_corr = 0.0, E_sd = 0.0, t_diss_on = 0.125, - e_kin_corr = 0.0, e_pot_corr = 0.0, e_pot_BH_corr = 0.0, - e_summ_corr = 0.0, - mcm, rcm_mod, vcm_mod, - rcm_sum=0.0, vcm_sum=0.0, - eps=0.0, eps2, - xcm[3], vcm[3], mom[3], - xdc[3], vdc[2], - over2=(1.0/2.0), over3=(1.0/3.0), over6=(1.0/6.0), - a2_mod, adot2_mod, - dt_tmp, dt2half, dt3over6, dt4over24, dt5over120, - dtinv, dt2inv, dt3inv, - a0mia1, ad04plad12, ad0plad1, - a2[3], a3[3], a2dot1[3], - a1abs, adot1abs, a2dot1abs, a3dot1abs, - Timesteps=0.0, n_act_sum=0.0, n_act_distr[N_MAX], g6_calls=0.0, g6_calls_sum=0.0, - tmp, tmp_r, tmp_v, tmp_rv, tmp_cpu, - tmp_pot, tmp_a, tmp_adot, - tmp_a_bh, tmp_adot_bh, - tmp_a_bh_pn0, tmp_a_bh_pn1, tmp_a_bh_pn2, tmp_a_bh_pn2_5, tmp_a_bh_pn3, tmp_a_bh_pn3_5, tmp_a_bh_spin, - tmp_adot_bh_pn0, tmp_adot_bh_pn1, tmp_adot_bh_pn2, tmp_adot_bh_pn2_5, tmp_adot_bh_pn3, tmp_adot_bh_pn3_5, tmp_adot_bh_spin; - -double R[3], V[3], L[3], mju, dr2, dr, dv2, dv, dl2, dl, pot_eff; - -char processor_name[MPI_MAX_PROCESSOR_NAME], - inp_fname[30], inp_fname_gmc[30], - out_fname[30], dbg_fname[30]; - -/* global variables */ - -int N, N_star, N_gmc, N_bh, - ind[N_MAX], name[N_MAX]; -double m[N_MAX], x[N_MAX][3], v[N_MAX][3], - pot[N_MAX], a[N_MAX][3], adot[N_MAX][3], - t[N_MAX], dt[N_MAX]; - -/* local variables */ - -int n_loc, ind_loc[N_MAX_loc]; -double m_loc[N_MAX_loc], x_loc[N_MAX_loc][3], v_loc[N_MAX_loc][3], - pot_loc[N_MAX_loc], a_loc[N_MAX_loc][3], adot_loc[N_MAX_loc][3], - t_loc[N_MAX_loc], dt_loc[N_MAX_loc]; - -/* data for active particles */ - -int n_act, ind_act[N_MAX]; -double m_act[N_MAX], - x_act[N_MAX][3], v_act[N_MAX][3], - pot_act[N_MAX], a_act[N_MAX][3], adot_act[N_MAX][3], - t_act[N_MAX], dt_act[N_MAX], - x_act_new[N_MAX][3], v_act_new[N_MAX][3], - pot_act_new[N_MAX], a_act_new[N_MAX][3], adot_act_new[N_MAX][3], - pot_act_tmp[N_MAX], a_act_tmp[N_MAX][3], adot_act_tmp[N_MAX][3], - pot_act_tmp_loc[N_MAX], a_act_tmp_loc[N_MAX][3], adot_act_tmp_loc[N_MAX][3]; - -FILE *inp, *out, *tmp_file, *dbg; - -double CPU_time_real0, CPU_time_user0, CPU_time_syst0; -double CPU_time_real, CPU_time_user, CPU_time_syst; - - -#ifdef TIMING - -double CPU_tmp_real0, CPU_tmp_user0, CPU_tmp_syst0; -double CPU_tmp_real, CPU_tmp_user, CPU_tmp_syst; - -double DT_TOT, - DT_ACT_DEF1, DT_ACT_DEF2, DT_ACT_DEF3, DT_ACT_PRED, - DT_ACT_GRAV, DT_EXT_GRAV, - DT_GMC_GRAV, DT_GMC_GMC_GRAV, DT_EXT_GMC_GRAV, - DT_ACT_CORR, DT_ACT_LOAD, - DT_STEVOL, DT_STARDISK, DT_STARDESTR; - -double DT_ACT_REDUCE; - -#endif - -/* some local settings for G6a board's */ - -int clusterid, ii, nn, numGPU; - -//#ifdef SAPPORO -#if defined(SAPPORO) || defined(GUINNESS) - -int npipe=GPU_PIPE, index_i[GPU_PIPE]; -double h2_i[GPU_PIPE], x_i[GPU_PIPE][3], v_i[GPU_PIPE][3], - p_i[GPU_PIPE], a_i[GPU_PIPE][3], jerk_i[GPU_PIPE][3]; -double new_tunit=51.0, new_xunit=51.0; - -#elif defined YEBISU - -int npipe=G6_NPIPE, index_i[G6_NPIPE]; -double h2_i[G6_NPIPE], x_i[G6_NPIPE][3], v_i[G6_NPIPE][3], - p_i[G6_NPIPE], a_i[G6_NPIPE][3], jerk_i[G6_NPIPE][3]; -int new_tunit=51, new_xunit=51; -#else - -int npipe=48, index_i[48]; -double h2_i[48], x_i[48][3], v_i[48][3], - p_i[48], a_i[48][3], jerk_i[48][3]; -int new_tunit=51, new_xunit=51; - -#endif - -int aflag=1, jflag=1, pflag=1; - -double ti=0.0, a2by18[3], a1by6[3], aby2[3]; - -/* normalization... */ - -#ifdef NORM -double m_norm, r_norm, v_norm, t_norm; -#endif - -double eps_BH=0.0; - -/* external potential... */ - -#ifdef EXTPOT - -#ifdef EXTPOT_GAL -double m_bulge, a_bulge, b_bulge, - m_disk, a_disk, b_disk, - m_halo, a_halo, b_halo, -// x_ext, y_ext, z_ext, -// vx_ext, vy_ext, vz_ext, -// x_ij, y_ij, z_ij, -// vx_ij, vy_ij, vz_ij, rv_ij, - x2_ij, y2_ij, z2_ij, - r_tmp, r2_tmp, z_tmp, z2_tmp; -#endif - -#ifdef EXTPOT_GAL_DEH -double m_ext, r_ext, g_ext, - tmp_r2, tmp_r3, dum, dum2, dum3, dum_g, tmp_g; -#endif - -#ifdef EXTPOT_GAL_LOG -double v_halo, r_halo, - v2_halo, r2_halo, r2_r2_halo, - x2_ij, y2_ij, z2_ij; -#endif - -#ifdef EXTPOT_BH -double r2, rv_ij, - m_bh, b_bh, eps_bh; -#endif - -#ifdef EXTPOT_SC -int M_SC=M_SC_DIM; -double r_sc[M_SC_DIM], m_sc[M_SC_DIM], p_sc[M_SC_DIM]; -double M_R, pot_out_R; -#endif - - -double pot_ext[N_MAX], pot_act_ext[N_MAX]; // for all EXTPOT - -#endif - - -int i_bh, i_bh1, i_bh2, - num_bh = 0, num_bh1 = 0, num_bh2 = 0; - -double m_bh, m_bh1, m_bh2, b_bh, - r, r2, - x_ij, y_ij, z_ij, - vx_ij, vy_ij, vz_ij, rv_ij; - -#ifdef STEVOL -int num_dead, event[N_MAX]; -double dt_stevol, t_stevol; -double m0[N_MAX], ZZZ[N_MAX], t0[N_MAX]; -#endif - - -#ifdef STEVOL_SSE - -int num_dead, event[N_MAX], event_old[N_MAX]; -double dt_stevol, t_stevol; -double m0[N_MAX], ZZZ[N_MAX], t0[N_MAX]; -double SSEMass[N_MAX], SSESpin[N_MAX], SSERad[N_MAX], - SSELum[N_MAX], SSETemp[N_MAX]; -double SSEMV[N_MAX], SSEBV[N_MAX]; - -double Rgal = 10000.0; //Distance of star from sun for artificial CMD with observational errors [pc] - -int van_kick=0; - -extern void zcnsts_(double *z, double *zpars); -extern void evolv1_(int *kw, double *mass, double *mt, double *r, double *lum, double *mc, double *rc, double *menv, double *renv, double *ospin, double *epoch, double *tms, double *tphys, double *tphysf, double *dtp, double *z, double *zpars, double *vkick, double *vs); -extern int cmd(int kstar, double rstar, double lstar, double Rgal, double *abvmag, double *vmag, double *BV, double *Teff, double *dvmag, double *dBV); - -struct{ - double neta; - double bwind; - double hewind; - double mxns; -} value1_; - -struct{ - double pts1; - double pts2; - double pts3; -} points_; - -struct{ - double sigma; - int bhflag; -} value4_; - -struct{ - int idum; -} value3_; - -struct{ - int ceflag; - int tflag; - int ifflag; - int nsflag; - int wdflag; -} flags_; - -struct{ - double beta; - double xi; - double acc2; - double epsnov; - double eddfac; - double gamma; -} value5_; - -struct{ - double alpha1; - double lambda; -} value2_; - -#include "cmd.c" - -#endif - - - -#ifdef BBH_INF -int inf_event[N_MAX]; -double x_bbhc[3], v_bbhc[3], DR2, tmp_r2; -double DV2, EB, SEMI_a, SEMI_a2; -#endif - - - -/* STARDISK Drag force... */ - -#ifdef STARDISK - -double a_drag[N_MAX][3], adot_drag[N_MAX][3]; -//double a_drag[N_MAX][3], adot_drag[N_MAX][3]; - -double z_new_drag, r_new_drag2, hz; - -#ifdef SPH_GAS - -#define N_GAS_MAX (128*KB) -#define NB 50 - -int N_GAS, ind_gas[N_GAS_MAX], sosed_gas[N_GAS_MAX][NB]; -double m_gas[N_GAS_MAX], x_gas[N_GAS_MAX][3], v_gas[N_GAS_MAX][3], h_gas[N_GAS_MAX], DEN_gas[N_GAS_MAX]; - -//int sosed[NB]; -double d[N_MAX]; // podrazumevajem chto: N_MAX > N_GAS_MAX vsegda ! -double x_star, y_star, z_star, DEN_star; -double h_min = 1.0E-04, h_max = 1.0; - -#include "def_H.c" -#include "def_DEN.c" -#include "def_DEN_STAR.c" - -#endif // SPH_GAS - -#include "drag_force.c" - -// calculate the SG for the set of active particles which is deepley INSIDE the EXT BH infl. rad. -// EXPERIMENTAL feature !!! - -//int SG_CALC = 0; -//double summ_tmp_r = 0.0, aver_tmp_r = 0.0, R_INT_CUT = 1.0E-03; - -#endif // STARDISK - - - -/* GMC data... */ - -#ifdef GMC - -double dt_gmc, E_pot_gmc, E_pot_gmc_gmc, E_pot_ext_gmc, E_kin_gmc; - -#define N_gmc_MAX (500) - -int N_gmc, ind_gmc[N_gmc_MAX], name_gmc[N_gmc_MAX]; - -double m_gmc[N_gmc_MAX], pot_gmc_gmc[N_gmc_MAX], pot_ext_gmc[N_gmc_MAX], - x_gmc[N_gmc_MAX][3], v_gmc[N_gmc_MAX][3], a_gmc[N_gmc_MAX][3], - eps_gmc[N_gmc_MAX]; - -double pot_gmc[N_MAX]; - -#endif // GMC - - -/****************************************************************************/ - -/****************************************************************************/ -#ifdef ADD_N_BH - -double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3]; - -double pot_bh1, a_bh1[3], adot_bh1[3], - pot_bh2, a_bh2[3], adot_bh2[3]; - -//double eps_BH = 0.0; - -#include "n_bh.c" - -/* -int calc_force_n_BH(double m1, double xx1[], double vv1[], - double m2, double xx2[], double vv2[], - double eps_BH, - double pot_n1, double a_n1[], double adot_n1[], - double pot_n2, double a_n2[], double adot_n2[]) -*/ - -/* - INPUT - -m1 - mass of the 1 BH -xx1[0,1,2] - coordinate of the 1 BH -vv1[0,1,2] - velocity of the 1 BH - -m2 - mass of the 2 BH -xx2[0,1,2] - coordinate of the 2 BH -vv2[0,1,2] - velocity of the 2 BH - -eps_BH - force softening, can be even exactly 0.0 ! - - OUTPUT - -pot_n1 for the 1 BH -a_n1 [0,1,2] for the 1 BH -adot_n1 [0,1,2] for the 1 BH - -pot_n2 for the 2 BH -a_n2 [0,1,2] for the 2 BH -adot_n2 [0,1,2] for the 2 BH - -return - 0 if everything OK -*/ - -#endif -/****************************************************************************/ - -/****************************************************************************/ -#ifdef ADD_PN_BH - -double C_NB = 477.12; - -int usedOrNot[7] = {1, 1, 1, 1, 0, 0, 0}; - -double a_pn1[7][3], adot_pn1[7][3], - a_pn2[7][3], adot_pn2[7][3]; - -double s_bh1[3] = {0.0, 0.0, 1.0}; -double s_bh2[3] = {0.0, 0.0, 1.0}; - -#ifdef ADD_SPIN_BH // eto rabotajet vsegda !!! -#include "pn_bh_spin.c" -#else -#include "pn_bh.c" // eto staraja versija, bolshe ne rabotajet !!! -#endif - -/* -int calc_force_pn_BH(double m1, double xx1[], double vv1[], double ss1[], - double m2, double xx2[], double vv2[], double ss2[], - double CCC_NB, double dt_bh, - int usedOrNot[], - double a_pn1[][3], double adot_pn1[][3], - double a_pn2[][3], double adot_pn2[][3]) -*/ - -/* - INPUT - -m1 - mass of the 1 BH -xx1[0,1,2] - coordinate of the 1 BH -vv1[0,1,2] - velocity of the 1 BH -spin1[0,1,2] - normalized spin of the 1 BH - -m2 - mass of the 2 BH -xx2[0,1,2] - coordinate of the 2 BH -vv2[0,1,2] - velocity of the 2 BH -spin2[0,1,2] - normalized spin of the 2 BH - -CCC_NB - Speed of light "c" in internal units -dt_BH - timestep of the BH's, needed for the SPIN integration - -usedOrNot[PN0, PN1, PN2, PN2.5, PN3, PN3.5, SPIN] - different PN term usage: PN1, PN2, PN2.5, PN3, PN3.5, SPIN - 0 1 2 3 4 5 6 - - OUTPUT - -a_pn1 [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 1 BH -adot_pn1[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 1 BH - -a_pn2 [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 2 BH -adot_pn2[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 2 BH - -return - 0 if everything OK - - 505 if BH's separation < 4 x (RSwarch1 + RSwarch2) -*/ - -#endif - -/****************************************************************************/ - -/****************************************************************************/ -/* RAND_MAX = 2147483647 */ -/* my_rand : 0.0 - 1.0 */ -/* my_rand2 : -1.0 - 1.0 */ - -double my_rand(void) -{ -return( (double)(rand()/(double)RAND_MAX) ); -} - -double my_rand2(void) -{ -return (double)(2.0)*((rand() - RAND_MAX/2)/(double)RAND_MAX); -} -/****************************************************************************/ - -#ifdef STARDESTR -#include "star_destr.c" -#endif - -#ifdef STARDESTR_EXT -#include "star_destr_ext.c" -#endif - -#ifdef ACT_DEF_LL -#include "act_def_linklist.c" -#endif - -#ifdef ETICS -double t_exp, dt_exp=ETICS_DTSCF; // t_exp is just the initial value -#ifdef ETICS_CEP -int grapite_cep_index; -#endif -#endif - -/****************************************************************************/ -void get_CPU_time(double *time_real, double *time_user, double *time_syst) - { - struct rusage xxx; - double sec_u, microsec_u, sec_s, microsec_s; - struct timeval tv; - - - getrusage(RUSAGE_SELF,&xxx); - - sec_u = xxx.ru_utime.tv_sec; - sec_s = xxx.ru_stime.tv_sec; - - microsec_u = xxx.ru_utime.tv_usec; - microsec_s = xxx.ru_stime.tv_usec; - - *time_user = sec_u + microsec_u * 1.0E-06; - *time_syst = sec_s + microsec_s * 1.0E-06; - -// *time_real = time(NULL); - - gettimeofday(&tv, NULL); - *time_real = tv.tv_sec + 1.0E-06 * tv.tv_usec; - - *time_user = *time_real; - - } -/****************************************************************************/ - -#ifdef STEVOL -/****************************************************************************/ -double t_dead(double m, double Z) -/* -m - mass of the star (in Msol) -Z - metallicity (absolut values - i.e. Zsol = 0.0122) -t - return value of star lifetime (in Year) - -Raiteri, Villata & Navarro, 1996, A&A, 315, 105 -*/ -{ - double a0, a1, a2, lZ, lm, tmp=0.0; - - m *= (m_norm/Msol); - - lZ = log10(Z); - lm = log10(m); - - a0 = 10.130 + 0.07547*lZ - 0.008084*lZ*lZ; - a1 = -4.424 - 0.79390*lZ - 0.118700*lZ*lZ; - a2 = 1.262 + 0.33850*lZ + 0.054170*lZ*lZ; - - tmp = a0 + a1*lm + a2*lm*lm; - - tmp = pow(10.0,tmp); - - tmp *= (Year/t_norm); - - return(tmp); -} -/****************************************************************************/ - -/****************************************************************************/ -double m_loss_old(double m, double Z) -/* -m - initial mass of the star (in Msol) -Z - metallicity (e.g. Zsol = 0.02) -m_ret - returned maass to ISM from the star (in Msol) - -approx. from data of van den Hoek & Groenewegen, 1997, A&AS, 123, 305 -*/ -{ - double tmp=0.0; - - m *= (m_norm/Msol); - -// tmp = (-0.46168 + 0.912274 * m); // first approx. - - tmp = -0.420542*pow(Z, -0.0176601) + (0.901495 + 0.629441*Z) * m; - - tmp *= (Msol/m_norm); - - return(tmp); -} -/****************************************************************************/ - -/****************************************************************************/ -double m_curr_old(double m0, double Z, double t0, double t) -{ - double td, ml, tmp=0.0; - - td = t_dead(m0, Z); - ml = m_loss_old(m0, Z); - - if( (t-t0) < td ) - { - tmp = m0; // instant mass loss -// tmp = m0 - (t-t0)/(td-t0) * ml; // linear mass loss - } - else - { - tmp = m0 - ml; - - if(event[i] == 0) - { - num_dead++; - event[i] = 1; - } - - } /* if( (t-t0) < td ) */ - - return(tmp); -} -/****************************************************************************/ - -/****************************************************************************/ -double m_loss(double m, double Z) -/* -m - initial mass of the star (in norm. units) -Z - metallicity (absolut value - i.e. Zsol = 0.0122) -m_loss - returned maass to ISM from the star (in norm. units) -*/ -{ - double tmp=0.0; - - m *= (m_norm/Msol); - - if( m < 8.0) - { -// approx. from data of van den Hoek & Groenewegen, 1997, A&AS, 123, 305 - -// tmp = (-0.46168 + 0.912274 * m); // first approx. - tmp = -0.420542*pow(Z, -0.0176601) + (0.901495 + 0.629441*Z) * m; // second approx. - } - else - { -// approx. from data of Woosley & Weaver, 1995, ApJS, 110, 181 - -// tmp = 0.813956*pow(m, 1.04214); - tmp = 0.813956*pow(m, 1.04); - } - - tmp *= (Msol/m_norm); - - return(tmp); -} -/****************************************************************************/ - -/****************************************************************************/ -double m_curr(double m0, double Z, double t0, double t) -{ - double td, ml, tmp=0.0; - - tmp = m0; - - if( event[i] == 0 ) - { - - td = t_dead(m0, Z); - - if( (t-t0) > td ) - { - - ml = m_loss(m0, Z); - - tmp = m0 - ml; - - num_dead++; - -// wd http://en.wikipedia.org/wiki/Chandrasekhar_limit 1.38 - - if( tmp * (m_norm/Msol) < 1.38 ) - { - event[i] = 1; - } - else - { - -// ns http://en.wikipedia.org/wiki/Neutron_stars 1.38 - ~2.0 ??? - - if( tmp * (m_norm/Msol) > 1.38 ) event[i] = 2; - -// bh http://en.wikipedia.org/wiki/Tolman-Oppenheimer-Volkoff_limit ~1.5 - ~3.0 ??? - - if( tmp * (m_norm/Msol) > 2.00 ) event[i] = 3; - - } - - } /* if( (t-t0) > td ) */ - - } /* if( event[i] == 0 ) */ - - return(tmp); -} -/****************************************************************************/ -#endif - - -#ifdef STEVOL_SSE -/****************************************************************************/ -double m_curr(int ipart, double m0, double Z, double t0, double t) -{ - //set up parameters and variables for SSE (Hurley, Pols & Tout 2002) - int kw=0; //stellar type - double mass; //initial mass - double mt=0.0; //actual mass - double r=0.0; //radius - double lum=0.0; //luminosity - double mc=0.0; //core mass - double rc=0.0; //core radius - double menv=0.0; //envelope mass - double renv=0.0; //envelope radius - double ospin=0.0; //spin - double epoch=0.0; //time spent in current evolutionary state - double tms=0.0; //main-sequence lifetime - double tphys; //initial age - double tphysf; //final age - double dtp=0.0; //data store value, if dtp>tphys no data will be stored - double z; //metallicity - double zpars[20]; //metallicity parameters - double vkick=0.0; //kick velocity for compact remnants - double vs[3]; //kick velocity componets - - // M_V_[mag] V_[mag] B-V_[mag] T_eff_[K] dV_[mag] d(B-V) - double abvmag, vmag, BV, Teff, dvmag, dBV; - - mass = m0*(m_norm/Msol); // m0[i] initial mass of the stars Msol - mt = mass; // current mass for output in Msol - tphys = t0*(t_norm/Myr); // t0[i] time of the star formation Myr - tphysf = t*(t_norm/Myr); // t[i] current time for the star Myr - dtp = 0.0; // output parameter, just set to 0.0 !!! - z = Z; // ZZZ[i] of the star - -/* - if(ipart > 9995) - { - printf("Initial parameters: \n"); - printf("M = %g [Mo] \n", mass); - printf("Z = %g \n", z); - printf("kw = %d \n", kw); - printf("spin = %g \n", ospin); - printf("epoch = %g [Myr] \n", epoch); - printf("t_beg = %g [Myr] \n", tphys); - printf("t_end = %g [Myr] \n", tphysf); - printf("R = %g [R_o] \n", r); - printf("L = %g [L_o] \n", lum); - printf("V_kick = %g [km/s] \n", vkick); - } -*/ - - for (k=0; k<20; k++) zpars[k] = 0; - zcnsts_(&z,zpars); //get metallicity parameters - - evolv1_(&kw, &mass, &mt, &r, &lum, &mc, &rc, &menv, &renv, &ospin, &epoch, &tms, &tphys, &tphysf, &dtp, &z, zpars, &vkick, vs); - cmd(kw, r, lum, Rgal, &abvmag, &vmag, &BV, &Teff, &dvmag, &dBV); - -/* - if(ipart > 9995) - { - printf("Final parameters: \n"); - printf("M = %g [Mo] \n", mt); - printf("Z = %g \n", z); - printf("kw = %d \n", kw); - printf("spin = %g \n", ospin); - printf("epoch = %g [Myr] \n", epoch); - printf("t_beg = %g [Myr] \n", tphys); - printf("t_end = %g [Myr] \n", tphysf); - printf("R = %g [R_o] \n", r); - printf("L = %g [L_o] \n", lum); - printf("V_kick = %g [km/s] \n", vkick); - printf("M_V_[mag]\tV_[mag]\t\tB-V_[mag]\tT_eff_[K]\tdV_[mag]\td(B-V)\n"); - printf("%.3E\t%.3E\t%.3E\t%.3E\t%.3E\t%.3E\n", abvmag, vmag, BV, Teff, dvmag, dBV); -// printf("%g\t%g\t%g\t%g\t%g\t%g\n", abvmag, vmag, BV, Teff, dvmag, dBV); - } -*/ - - event[ipart] = kw; // Evolution type (0 - 15). - SSEMass[ipart] = mt; // Mass in Msol - SSESpin[ipart] = ospin; // Spin in SI ??? [kg*m*m/s] - SSERad[ipart] = r; // Radius in Rsol - SSELum[ipart] = lum; // Luminosity in Lsol - SSETemp[ipart] = Teff; // Temperature - - SSEMV[ipart] = abvmag; // M_V absolute V magnitude - SSEBV[ipart] = BV; // B-V color index - -/* -c STELLAR TYPES - KW -c -c 0 - deeply or fully convective low mass MS star -c 1 - Main Sequence star -c 2 - Hertzsprung Gap -c 3 - First Giant Branch -c 4 - Core Helium Burning -c 5 - First Asymptotic Giant Branch -c 6 - Second Asymptotic Giant Branch -c 7 - Main Sequence Naked Helium star -c 8 - Hertzsprung Gap Naked Helium star -c 9 - Giant Branch Naked Helium star -c 10 - Helium White Dwarf -c 11 - Carbon/Oxygen White Dwarf -c 12 - Oxygen/Neon White Dwarf -c 13 - Neutron Star -c 14 - Black Hole -c 15 - Massless Supernova -*/ - - if( van_kick == 1 ) // we have a kick ??? - { - - if( (event_old[ipart] < 10) && ( (event[ipart] == 13) || (event[ipart] == 14) ) ) // NS or BH - { - - if(myRank == rootRank) - { -// printf("KICK: %06d %.6E [Myr] %02d (%02d) %.6E [Mo] %.6E [km/s] [% .3E, % .3E, % .3E] OLD: [% .3E, % .3E, % .3E] %.6E [Mo] \n", ipart, tphysf, kw, event_old[ipart], mt, vkick, vs[0], vs[1], vs[2], v[ipart][0]*(v_norm/km), v[ipart][1]*(v_norm/km), v[ipart][2]*(v_norm/km), mass); - printf("KICK: %06d %.6E %02d %02d %.6E %.6E % .3E % .3E % .3E % .3E % .3E % .3E %.6E \n", ipart, tphysf, kw, event_old[ipart], mt, vkick, vs[0], vs[1], vs[2], v[ipart][0]*(v_norm/km), v[ipart][1]*(v_norm/km), v[ipart][2]*(v_norm/km), mass); - fflush(stdout); - } /* if(myRank == rootRank) */ - - v[ipart][0] += vs[0]*(km/v_norm); - v[ipart][1] += vs[1]*(km/v_norm); - v[ipart][2] += vs[2]*(km/v_norm); - } - - } // we have a kick ??? - - - event_old[ipart] = event[ipart]; - - - tmp = mt*(Msol/m_norm); // return the current mass(t-t0) of the star in NB units - - return(tmp); -} -/****************************************************************************/ -#endif - - -/****************************************************************************/ -void read_data() -{ - inp = fopen(inp_fname,"r"); - fscanf(inp,"%d \n", &diskstep); - -//#ifdef STEVOL -// fscanf(inp,"%d \n", &N); -//#endif -// fscanf(inp,"%d %d %d %d \n", &N, &N_bh, &N_gmc, &N_star); - - fscanf(inp,"%d \n", &N); - - fscanf(inp,"%lE \n", &time_cur); - - for(i=0; i= (mass_frac[k]*m_tot) ) - { - lagr_rad[k] = var_sort[j]; - k++; - } -#endif - } - - fclose(out); - -// write out lagr-rad.dat - -#ifdef LAGR_RAD - out = fopen("lagr-rad.dat","a"); - fprintf(out,"%.6E ", time_cur); - for(k=0; k 100.0) jerk_i[ii][0] = 100.0; - if(ABS(jerk_i[ii][1]) > 100.0) jerk_i[ii][1] = 100.0; - if(ABS(jerk_i[ii][2]) > 100.0) jerk_i[ii][2] = 100.0; - } - -#ifdef SAPPORO - g6calc_firsthalf_(&clusterid, &n_loc, &nn, index_i, x_i, v_i , a_i, jerk_i, p_i, &eps2, h2_i); - g6calc_lasthalf_(&clusterid, &n_loc, &nn, index_i, x_i, v_i, &eps2, h2_i, a_i, jerk_i, p_i); -#else - g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i); - g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i); -#endif - - for(ii=0; ii r) ) - { - r1 = r_sc[ii]; r2 = r_sc[ii+1]; - m1 = m_sc[ii]; m2 = m_sc[ii+1]; - p1 = p_sc[ii]; p2 = m_sc[ii+1]; - break; - } - } - - m_tmp = m1 + (r-r1)*(m2-m1)/(r2-r1); - p_tmp = p1 + (r-r1)*(p2-p1)/(r2-r1); - -// if(myRank == rootRank) -// { -// printf("\t \t \t %.6E %.6E \n", r1, r2); -// printf("\t \t \t %.6E %.6E \n", m1, m2); -// printf("\t \t \t %.6E %.6E \n", p1, p2); -// printf("\t \t \t %06d %.6E %.6E %.6E \n", ii, r, m_tmp, p_tmp); -// fflush(stdout); -// } /* if(myRank == rootRank) */ - - *M_R = m_tmp; - *pot_out_R = p_tmp; - -//exit(-1); - - } -#endif -/****************************************************************************/ - -/****************************************************************************/ -void calc_ext_grav_zero() -{ - -#ifdef EXTPOT - - /* Define the external potential for all particles on all nodes */ - - ni = n_act; - - for(i=0; i 100.0) jerk_i[ii][0] = 100.0; - if(ABS(jerk_i[ii][1]) > 100.0) jerk_i[ii][1] = 100.0; - if(ABS(jerk_i[ii][2]) > 100.0) jerk_i[ii][2] = 100.0; - } - - g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i); - g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i); - - g6_calls++; -*/ - - for(ii=0; ii R_LIMITS_ALL) && (m[i] != 0.0) ) - if( (tmp_r > R_LIMITS_ALL) && (tmp_rv > 0.0) ) - { - - for(k=0;k<3;k++) - { -// x[i][k] = R_LIMITS_ALL + 1.0*my_rand2(); - v[i][k] *= 1.0E-01; - } /* k */ - - } /* if */ - - } /* i */ - -#endif - - -#ifdef CMCORR - - /* possible CM correction of the initial datafile... */ - - if( (diskstep == 0) && (time_cur == 0.0) ) cm_corr(); -// if( (time_cur = 0.0) ) cm_corr(); - -#endif - - - printf("\n"); - printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc); - printf("\n"); -#ifdef GMC - printf("N = %06d \t N_GMC = %06d \t eps = %.6E \n", N, N_gmc, eps); -#else - printf("N = %07d \t eps = %.6E \n", N, eps); -#endif - printf("t_beg = %.6E \t t_end = %.6E \n", time_cur, t_end); - printf("dt_disk = %.6E \t dt_contr = %.6E \n", dt_disk, dt_contr); - printf("dt_bh = %.6E \n", dt_bh); - printf("eta = %.6E \n", eta); - printf("\n"); - -#ifdef NORM - printf("Normalization: \n"); - printf("\n"); - printf("m_norm = %.4E [Msol] r_norm = %.4E [pc] \n", m_norm/Msol, r_norm/pc); - printf("v_morm = %.4E [km/s] t_morm = %.4E [Myr] \n", v_norm/km, t_norm/Myr); - printf("\n"); -#endif - -#ifdef EXTPOT - printf("External Potential: \n"); - printf("\n"); - -#ifdef EXTPOT_GAL - printf("m_bulge = %.4E [Msol] a_bulge = %.4E b_bulge = %.4E [kpc] \n", m_bulge*m_norm/Msol, a_bulge*r_norm/kpc, b_bulge*r_norm/kpc); - printf("m_disk = %.4E [Msol] a_disk = %.4E b_disk = %.4E [kpc] \n", m_disk*m_norm/Msol, a_disk*r_norm/kpc, b_disk*r_norm/kpc); - printf("m_halo = %.4E [Msol] a_halo = %.4E b_halo = %.4E [kpc] \n", m_halo*m_norm/Msol, a_halo*r_norm/kpc, b_halo*r_norm/kpc); -#endif - -#ifdef EXTPOT_BH - printf("m_bh = %.4E b_bh = %.4E \n", m_bh, b_bh); -// printf("m_bh = %.4E [Msol] b_bh = %.4E [kpc] \n", m_bh*m_norm/Msol, b_bh*r_norm/kpc); -#endif - -#ifdef EXTPOT_GAL_DEH - printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n", m_ext, r_ext, g_ext); -#endif - -#ifdef EXTPOT_GAL_LOG - printf("v_halo = %.6E r_halo = %.6E \n", v_halo, r_halo); -#endif - -#ifdef EXTPOT_SC - read_SC_mass(); - - printf("EXTPOT_SC: # of points %06d \t M_SC_TOT = %.6E \n", M_SC_DIM-1, m_sc[M_SC_DIM-1]); -// printf("EXTPOT_SC: %.6E \t %.6E \n", m_sc[0], m_sc[M_SC_DIM-1]); -// printf("EXTPOT_SC: %.6E \t %.6E \n", p_sc[0], p_sc[M_SC_DIM-1]); -#endif - - printf("\n"); -#endif - - - - fflush(stdout); - - - if( (diskstep == 0) && (time_cur == 0.0) ) -// if( (time_cur = 0.0) ) - { - -#ifdef SPH_GAS - write_data_SPH(); -#endif - -#ifdef GMC - write_snap_GMC(); - write_cont_GMC(); -#endif - -// write_snap_data(); -// write_cont_data(); - - out = fopen("contr.dat","w"); - fclose(out); - -#ifdef TIMING - out = fopen("timing.dat","w"); - fclose(out); -#endif - -#ifdef ADD_BH2 - -#ifdef BH_OUT - out = fopen("bh.dat","w"); - fclose(out); -#endif - -#ifdef BH_OUT_NB - out = fopen("bh_nb.dat","w"); - fclose(out); -#endif - -#else - -#ifdef ADD_BH1 - -#ifdef BH_OUT - out = fopen("bh.dat","w"); - fclose(out); -#endif - -#ifdef BH_OUT_NB - out = fopen("bh_nb.dat","w"); - fclose(out); -#endif - -#endif // ADD_BH1 - -#endif // ADD_BH2 - - -#ifdef STARDESTR - -#ifdef ADD_BH2 - - out = fopen("mass-bh.dat","w"); - m_bh1 = m[0]; - m_bh2 = m[1]; - num_bh1 = 0; - num_bh2 = 0; - fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", - time_cur, m_bh1, num_bh1, m_bh2, num_bh2); - fclose(out); - - out = fopen("mass-bh.con","w"); - m_bh1 = m[0]; - m_bh2 = m[1]; - num_bh1 = 0; - num_bh2 = 0; - fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", - time_cur, m_bh1, num_bh1, m_bh2, num_bh2); - fclose(out); - -#else - -#ifdef ADD_BH1 - - out = fopen("mass-bh.dat","w"); - m_bh1 = m[0]; - num_bh1 = 0; - fprintf(out,"%.8E \t %.8E %06d \n", - time_cur, m_bh1, num_bh1); - fclose(out); - - out = fopen("mass-bh.con","w"); - m_bh1 = m[0]; - num_bh1 = 0; - fprintf(out,"%.8E \t %.8E %06d \n", - time_cur, m_bh1, num_bh1); - fclose(out); - -#endif // ADD_BH1 - -#endif // ADD_BH2 - -#endif // STARDESTR - - -#ifdef STARDESTR_EXT - - num_bh = 0; - out = fopen("mass-bh.dat","w"); - fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); - fclose(out); - - num_bh = 0; - out = fopen("mass-bh.con","w"); - fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); - fclose(out); - -#ifdef BH_OUT_NB_EXT - out = fopen("bh_nb.dat","w"); - fclose(out); -#endif - -#endif // STARDESTR_EXT - - -#ifdef DEBUG_extra -#ifdef LAGR_RAD - out = fopen("lagr-rad.dat","w"); - fclose(out); -#endif -#endif - - -#ifdef BBH_INF - out = fopen("bbh.inf","w"); - fclose(out); -#endif - - - } - else // if(diskstep == 0) - { - -//#ifdef LAGR_RAD -// out = fopen("lagr-rad.dat","a"); -// fclose(out); -//#endif - -/* -#ifdef SPH_GAS - write_data_SPH(); -#endif - -#ifdef GMC - write_snap_GMC(); - write_cont_GMC(); -#endif - - write_snap_data(); - write_cont_data(); -*/ - -#ifdef STARDESTR - -#ifdef ADD_BH2 - -/* - inp = fopen("mass-bh.dat","r"); - while(feof(inp) == NULL) - { - fscanf(inp,"%lE %lE %d %lE %d", &tmp, &m_bh1, &num_bh1, &m_bh2, &num_bh2); - } - fclose(inp); - printf("%.8E \t %.8E %06d \t %.8E %06d \n", tmp, m_bh1, num_bh1, m_bh2, num_bh2); - fflush(stdout); -*/ - - inp = fopen("mass-bh.con","r"); - fscanf(inp,"%lE %lE %d %lE %d", &tmp, &m_bh1, &num_bh1, &m_bh2, &num_bh2); - fclose(inp); - printf("%.8E \t %.8E %06d \t %.8E %06d \n", tmp, m_bh1, num_bh1, m_bh2, num_bh2); - fflush(stdout); - -#else - -#ifdef ADD_BH1 - -/* - inp = fopen("mass-bh.dat","r"); - while(feof(inp) == NULL) - { - fscanf(inp,"%lE %lE %d", &tmp, &m_bh1, &num_bh1); - } - fclose(inp); - printf("%.8E \t %.8E \t %06d \n", tmp, m_bh1, num_bh1); - fflush(stdout); -*/ - - inp = fopen("mass-bh.con","r"); - fscanf(inp,"%lE %lE %d", &tmp, &m_bh1, &num_bh1); - fclose(inp); - printf("%.8E \t %.8E \t %06d \n", tmp, m_bh1, num_bh1); - fflush(stdout); - -#endif // ADD_BH1 - -#endif // ADD_BH2 - -#endif // STARDESTR - - -#ifdef STARDESTR_EXT - -/* - inp = fopen("mass-bh.dat","r"); - while(feof(inp) == NULL) - { - fscanf(inp,"%lE %lE %d", &tmp, &m_bh, &num_bh); - } - fclose(inp); - printf("%.8E \t %.8E \t %06d \n", tmp, m_bh, num_bh); - fflush(stdout); -*/ - - inp = fopen("mass-bh.con","r"); - fscanf(inp,"%lE %lE %d", &tmp, &m_bh, &num_bh); - fclose(inp); - printf("%.8E \t %.8E \t %06d \n", tmp, m_bh, num_bh); - fflush(stdout); - -#endif // STARDESTR_EXT - - - } // if(diskstep == 0) - - - get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0); - - } /* if(myRank == rootRank) */ - - - - -#ifdef STEVOL_SSE -//SSE internal parameters (see Hurley, Pols & Tout 2000) - -value1_.neta = 0.5; //Reimers mass-loss coefficent (neta*4x10^-13; 0.5 normally) -value1_.bwind = 0.0; //Binary enhanced mass loss parameter (inactive for single) -value1_.hewind = 1.0; //Helium star mass loss factor (1.0 normally) -value1_.mxns = 3.0; //Maximum NS mass (1.8, nsflag=0; 3.0, nsflag=1) - -points_.pts1 = 0.05; //Time-step parameter in evolution phase: MS (0.05) -points_.pts2 = 0.01; //Time-step parameter in evolution phase: GB, CHeB, AGB, HeGB (0.01) -points_.pts3 = 0.02; //Time-step parameter in evolution phase: HG, HeMS (0.02) - -//value4_.sigma = 190.0; //Kick velocities - standard values... -value4_.sigma = 265.0; //Kick velocities - Hobbs++2005 values... -value4_.bhflag = 1; //bhflag > 0 allows velocity kick at BH formation - -value3_.idum = 19640916; // random number seed !!! - -//BSE internal parameters (see Hurley, Pols & Tout 2002) - -flags_.ceflag = 0; //ceflag > 0 activates spin-energy correction in common-envelope (0) #defunct# -flags_.tflag = 1; //tflag > 0 activates tidal circularisation (1) -flags_.ifflag = 0; //ifflag > 0 uses WD IFMR of HPE, 1995, MNRAS, 272, 800 (0) -flags_.nsflag = 1; //nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ, 572, 407 (1) -flags_.wdflag = 1; //wdflag > 0 uses modified-Mestel cooling for WDs (0) - -value5_.beta = 0.125; //beta is wind velocity factor: proportional to vwind**2 (1/8) -value5_.xi = 1.0; //xi is the wind accretion efficiency factor (1.0) -value5_.acc2 = 1.5; //acc2 is the Bondi-Hoyle wind accretion factor (3/2) -value5_.epsnov = 0.001; //epsnov is the fraction of accreted matter retained in nova eruption (0.001) -value5_.eddfac = 1.0; //eddfac is Eddington limit factor for mass transfer (1.0) -value5_.gamma = -1.0; //gamma is the angular momentum factor for mass lost during Roche (-1.0) - -value2_.alpha1 = 1.0; //alpha1 is the common-envelope efficiency parameter (1.0) -value2_.lambda = 0.5; //lambda is the binding energy factor for common envelope evolution (0.5) - -#endif // STEVOL_SSE - - - /* Wait to all processors to finish his works... */ - MPI_Barrier(MPI_COMM_WORLD); - - /* Broadcast all useful values to all processors... */ - MPI_Bcast(&N, 1, MPI_INT, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&eps, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&eta, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&t_end, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&dt_disk, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&dt_contr, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&dt_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&time_cur, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - -#ifdef NORM - MPI_Bcast(&m_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&r_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&v_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&t_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); -#endif // NORM - -#ifdef EXTPOT - -#ifdef EXTPOT_GAL - MPI_Bcast(&m_bulge, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&a_bulge, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&b_bulge, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - - MPI_Bcast(&m_disk, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&a_disk, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&b_disk, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - - MPI_Bcast(&m_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&a_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&b_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - -/* - MPI_Bcast(&x_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&y_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&z_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - - MPI_Bcast(&vx_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&vy_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&vz_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); -*/ -#endif - -#ifdef EXTPOT_BH - MPI_Bcast(&m_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&b_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); -#endif - -#ifdef EXTPOT_GAL_DEH - MPI_Bcast(&m_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&r_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&g_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); -#endif - -#ifdef EXTPOT_GAL_LOG - MPI_Bcast(&v_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&r_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - - MPI_Bcast(&v2_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&r2_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); -#endif - -#endif // EXTPOT - - -#ifdef STARDESTR - MPI_Bcast(&m_bh1, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&num_bh1, 1, MPI_INT, rootRank, MPI_COMM_WORLD); - - MPI_Bcast(&m_bh2, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&num_bh2, 1, MPI_INT, rootRank, MPI_COMM_WORLD); -#endif // STARDESTR - - -#ifdef STARDESTR_EXT - MPI_Bcast(&num_bh, 1, MPI_INT, rootRank, MPI_COMM_WORLD); -#endif // STARDESTR_EXT - - - /* Wait to all processors to finish his works... */ - MPI_Barrier(MPI_COMM_WORLD); - - - - eta_s = eta/ETA_S_CORR; - eta_bh = eta/ETA_BH_CORR; - - eps2 = SQR(eps); - - dt_min = 1.0*pow(2.0, DTMINPOWER); - dt_max = 1.0*pow(2.0, DTMAXPOWER); - - if(dt_disk == dt_contr) - dt_max = dt_contr; - else - dt_max = MIN(dt_disk, dt_contr); - - if(dt_max > 1.0) dt_max = 1.0; - - t_disk = dt_disk*(1.0+floor(time_cur/dt_disk)); - t_contr = dt_contr*(1.0+floor(time_cur/dt_contr)); - t_bh = dt_bh*(1.0+floor(time_cur/dt_bh)); - - -if(myRank == rootRank) - { - printf("t_disk = %.6E t_contr = %.6E t_bh = %.6E \n", t_disk, t_contr, t_bh); - printf("\n"); - fflush(stdout); - } /* if(myRank == rootRank) */ - -/* - t_disk = time_cur + dt_disk; - t_contr = time_cur + dt_contr; - t_bh = time_cur + dt_bh; -*/ - -#ifdef STEVOL - dt_stevol = dt_max; - t_stevol = dt_stevol*(1.0 + floorl(time_cur/dt_stevol)); - -if(myRank == rootRank) - { - printf("t_stevol = %.6E dt_stevol = %.6E \n", t_stevol, dt_stevol); - printf("\n"); - fflush(stdout); - } /* if(myRank == rootRank) */ -#endif - -#ifdef STEVOL_SSE - dt_stevol = dt_max; - t_stevol = dt_stevol*(1.0 + floorl(time_cur/dt_stevol)); - -if(myRank == rootRank) - { - printf("t_stevol = %.6E dt_stevol = %.6E \n", t_stevol, dt_stevol); - printf("\n"); - fflush(stdout); - } /* if(myRank == rootRank) */ -#endif - - - for(i=0; i 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[i_bh1] -= pot_bh1; - pot[i_bh2] -= pot_bh2; - - for(k=0;k<3;k++) - { - a[i_bh1][k] -= a_bh1[k]; - a[i_bh2][k] -= a_bh2[k]; - - adot[i_bh1][k] -= adot_bh1[k]; - adot[i_bh2][k] -= adot_bh2[k]; - } - -// 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[i_bh1] += pot_bh1; - pot[i_bh2] += pot_bh2; - - for(k=0;k<3;k++) - { - a[i_bh1][k] += a_bh1[k]; - a[i_bh2][k] += a_bh2[k]; - - adot[i_bh1][k] += adot_bh1[k]; - adot[i_bh2][k] += adot_bh2[k]; - } - -#endif // ADD_N_BH - - -#ifdef ADD_PN_BH - -// calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk - - -/* - printf("PN: % .8E \t % .8E \t % .8E \n", x_bh1[0], x_bh1[1], x_bh1[2]); - printf("PN: % .8E \t % .8E \t % .8E \n", x_bh2[0], x_bh2[1], x_bh2[2]); - fflush(stdout); -*/ - - dt_bh_tmp = dt[0]; - - tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1, - m_bh2, x_bh2, v_bh2, s_bh2, - C_NB, dt_bh_tmp, usedOrNot, - a_pn1, adot_pn1, - a_pn2, adot_pn2); - - for(k=0;k<3;k++) - { - a[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k]; - a[i_bh2][k] += a_pn2[1][k] + a_pn2[2][k] + a_pn2[3][k] + a_pn2[4][k] + a_pn2[5][k] + a_pn2[6][k]; - - adot[i_bh1][k] += adot_pn1[1][k] + adot_pn1[2][k] + adot_pn1[3][k] + adot_pn1[4][k] + adot_pn1[5][k] + adot_pn1[6][k]; - adot[i_bh2][k] += adot_pn2[1][k] + adot_pn2[2][k] + adot_pn2[3][k] + adot_pn2[4][k] + adot_pn2[5][k] + adot_pn2[6][k]; - } - - if(myRank == rootRank) - { - if(tmp_i == 505) - { - printf("PN RSDIST: %.8E \t %.8E \n", Timesteps, time_cur); - fflush(stdout); - exit(-1); - } - } - -#endif // ADD_PN_BH - -#endif // ADD_BH2 - - - -#if defined(EXTPOT) || defined(GMC) - calc_ext_grav_zero(); -#endif - -#ifdef GMC - calc_gmc_self_grav(); - calc_ext_gmc_grav(); -#endif - - - /* Wait to all processors to finish his works... */ - MPI_Barrier(MPI_COMM_WORLD); - - - -// sjuda stavim DRAG... - -#ifdef STARDISK - - for(i=0; i dt_max) dt_tmp = dt_max; - -// dt_tmp = dt_min; - - -#ifdef STARDESTR - if(m[i] == 0.0) dt_tmp = dt_max; -#endif - -#ifdef STARDESTR_EXT - if(m[i] == 0.0) dt_tmp = dt_max; -#endif - - dt[i] = dt_tmp; - - -#ifdef DT_MIN_WARNING - if(myRank == 0) - { - if(dt[i] == dt_min) - { - printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[i], ind[i]); - fflush(stdout); - } - } -#endif - - } /* i */ - - - - -#ifdef ADD_BH2 - -/* define the min. dt over all the part. and set it also for the BH... */ - - min_dt = dt[0]; - - for(i=1; i 0.1) dt[i] = min_dt; - } /* i */ - -#endif // GMC2 - - -#ifdef GMC2222 - -/* define the max. dt for the GMC... */ - - for(i=1; i 0.1) dt[i] = dt_max; - } /* i */ - -#endif // GMC2 - - - -#ifdef ACT_DEF_LL - CreateLinkList(); -#ifdef DEBUG111 - if( myRank == rootRank ) check_linklist(0); -#endif -#endif - - - - /* Wait to all processors to finish his works... */ - MPI_Barrier(MPI_COMM_WORLD); - - /* Scatter the "local" vectors from "global" */ - MPI_Scatter(dt, n_loc, MPI_DOUBLE, dt_loc, n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - - /* Wait to all processors to finish his works... */ - MPI_Barrier(MPI_COMM_WORLD); - - /* load the new values for particles to the local GRAPE's */ - - nj=n_loc; - - /* load the nj particles to the G6 */ - - for(j=0; j 0) - for(int i=0; i= t_contr) -// { - MPI_Allreduce(pot_act_tmp, pot_act_new, n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); -// } /* if(min_t >= t_contr) */ - -// MPI_Allreduce(pot_act_tmp, pot_act_new, n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - - MPI_Allreduce(a_act_tmp, a_act_new, 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(adot_act_tmp, adot_act_new, 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - - /* Wait to all processors to finish his works... */ - MPI_Barrier(MPI_COMM_WORLD); - -#ifdef TIMING - get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst); - DT_ACT_REDUCE += (CPU_tmp_user - CPU_tmp_user0); -#endif - - - -#ifdef ADD_BH2 - -#ifdef ADD_N_BH - - m_bh1 = m_act[i_bh1]; - m_bh2 = m_act[i_bh2]; - - for(k=0;k<3;k++) - { - x_bh1[k] = x_act_new[i_bh1][k]; - v_bh1[k] = v_act_new[i_bh1][k]; - - x_bh2[k] = x_act_new[i_bh2][k]; - v_bh2[k] = v_act_new[i_bh2][k]; - } - -// 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; - - for(k=0;k<3;k++) - { - a_act_new[i_bh1][k] -= a_bh1[k]; - a_act_new[i_bh2][k] -= a_bh2[k]; - - adot_act_new[i_bh1][k] -= adot_bh1[k]; - adot_act_new[i_bh2][k] -= adot_bh2[k]; - } - -// 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; - - for(k=0;k<3;k++) - { - a_act_new[i_bh1][k] += a_bh1[k]; - a_act_new[i_bh2][k] += a_bh2[k]; - - adot_act_new[i_bh1][k] += adot_bh1[k]; - adot_act_new[i_bh2][k] += adot_bh2[k]; - } - -#endif // ADD_N_BH - - -#ifdef ADD_PN_BH - -// calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk - - dt_bh_tmp = dt[0]; - - tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1, - m_bh2, x_bh2, v_bh2, s_bh2, - C_NB, dt_bh_tmp, usedOrNot, - a_pn1, adot_pn1, - a_pn2, adot_pn2); - - for(k=0;k<3;k++) - { - a_act_new[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k]; - a_act_new[i_bh2][k] += a_pn2[1][k] + a_pn2[2][k] + a_pn2[3][k] + a_pn2[4][k] + a_pn2[5][k] + a_pn2[6][k]; - - adot_act_new[i_bh1][k] += adot_pn1[1][k] + adot_pn1[2][k] + adot_pn1[3][k] + adot_pn1[4][k] + adot_pn1[5][k] + adot_pn1[6][k]; - adot_act_new[i_bh2][k] += adot_pn2[1][k] + adot_pn2[2][k] + adot_pn2[3][k] + adot_pn2[4][k] + adot_pn2[5][k] + adot_pn2[6][k]; - } - - if(myRank == rootRank) - { - if(tmp_i == 505) - { - printf("PN RSDIST: TS = %.8E \t t = %.8E \n", Timesteps, time_cur); - fflush(stdout); - exit(-1); - } - } - -#endif // ADD_PN_BH - -#endif // ADD_BH2 - -#ifdef EXTPOT - calc_ext_grav(); -#endif - -#ifdef GMC - calc_gmc_self_grav(); - calc_ext_gmc_grav(); -#endif - - -// sjuda stavim DRAG... - -#ifdef STARDISK - -#ifdef TIMING - get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); -#endif - - calc_drag_force(); - - /* E_sd tut opjaty pljusujetsja... */ - - for(i=0; i dt_min) ) - { - power = log(dt_new)/log(2.0) - 1; -// power = log(dt_new)/log(2.0); - - dt_tmp = pow(2.0, (double)power); - } - - -#ifdef STARDISK - - // accretion disk criteria -// sjuda stavim izmenenija DT za schot DRAG pri priblizhenii k ploskosti Z=0 ... - - z_new_drag = x_act_new[i][2] + dt_tmp * v_act_new[i][2]; - r_new_drag2 = SQR(x_act_new[i][0]) + SQR(x_act_new[i][1]); - - if(r_new_drag2 > (R_CRIT*R_CRIT)) - { - hz = HZ; - } - else - { - hz = HZ*(sqrt(r_new_drag2)/R_CRIT); - } - -// if (fabs(x_act_new[i][2]) > HZ*2.2) - if (fabs(x_act_new[i][2]) > hz*2.2) - { - - if (z_new_drag * x_act_new[i][2] < 0.0) - { - dt_tmp *= 0.125; - } - else - { -// if (fabs(z_new_drag) < HZ*1.8) - if (fabs(z_new_drag) < hz*1.8) - { - dt_tmp *= 0.5; - } - } - - } - -#endif - - - if( (dt_new > 2.0*dt_tmp) && - (fmod(min_t, 2.0*dt_tmp) == 0.0) && - (2.0*dt_tmp <= dt_max) ) - { - dt_tmp *= 2.0; - } - - dt_act[i] = dt_tmp; - - t_act[i] = min_t; - - pot_act[i] = pot_act_new[i]; - - for(k=0; k<3; k++) - { - x_act[i][k] = x_act_new[i][k]; - v_act[i][k] = v_act_new[i][k]; - a_act[i][k] = a_act_new[i][k]; - adot_act[i][k] = adot_act_new[i][k]; - } /* k */ - -/* - if(myRank == rootRank) - { - tmp_r = sqrt( SQR(x_act_new[i][0]) + SQR(x_act_new[i][1]) + SQR(x_act_new[i][2]) ); - - if(ind_act[i] == 4232) - { - printf("SOS: TS = %.8E \t i = %06d \t n_act = %06d \t R = %.8E \n", Timesteps, ind_act[i], n_act, tmp_r); - printf(" x = % .6E % .6E % .6E \n", x_act_new[i][0], x_act_new[i][1], x_act_new[i][2]); - printf(" v = % .6E % .6E % .6E \n", v_act_new[i][0], v_act_new[i][1], v_act_new[i][2]); - printf(" a = % .6E % .6E % .6E \n", a_act_new[i][0], a_act_new[i][1], a_act_new[i][2]); - printf(" adot = % .6E % .6E % .6E \n", adot_act_new[i][0], adot_act_new[i][1], adot_act_new[i][2]); - printf(" t = %.6E %.6E \n", t_act[i], dt_act[i]); - fflush(stdout); -// exit(-1); - } - } -*/ - - -#ifdef DT_MIN_WARNING - if(myRank == 0) - { - if(dt_act[i] == dt_min) - { - printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]); - fflush(stdout); - } - } -#endif - - - - } /* i */ - - - - -/* define the min. dt over all the act. part. and set it also for the BH... */ - -#ifdef ADD_BH2 - - min_dt = dt_act[0]; - - for(i=1; i 0.1) dt_act[i] = min_dt; - } /* i */ - - -#endif // GMC2 - - -#ifdef GMC2222 - -/* define the max. dt for the GMC... */ - - for(i=1; i 0.1) dt_act[i] = dt_max; - } /* i */ - -#endif // GMC2 - - - - - - -#ifdef BBH_INF - if(myRank == rootRank) - { - - out = fopen("bbh.inf","a"); - - m_bh1 = m_act[i_bh1]; - m_bh2 = m_act[i_bh2]; - - for(k=0;k<3;k++) - { - x_bh1[k] = x_act[i_bh1][k]; - x_bh2[k] = x_act[i_bh2][k]; - - v_bh1[k] = v_act[i_bh1][k]; - v_bh2[k] = v_act[i_bh2][k]; - } - - for(k=0;k<3;k++) - { - x_bbhc[k] = (m_bh1*x_bh1[k] + m_bh2*x_bh2[k])/(m_bh1 + m_bh2); - v_bbhc[k] = (m_bh1*v_bh1[k] + m_bh2*v_bh2[k])/(m_bh1 + m_bh2); - } - - DR2 = SQR(x_bh1[0] - x_bh2[0]) + SQR(x_bh1[1] - x_bh2[1]) + SQR(x_bh1[2] - x_bh2[2]); - DV2 = SQR(v_bh1[0] - v_bh2[0]) + SQR(v_bh1[1] - v_bh2[1]) + SQR(v_bh1[2] - v_bh2[2]); - -// mju = m_bh1*m_bh2/(m_bh1 + m_bh2); - - EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2; - - SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB; - SEMI_a2 = SQR(SEMI_a); - - -// printf("INF: %.6E %.16E %.6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E \n", -// Timesteps, time_cur, SEMI_a, -// sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2], -// m_bh1, x_bh1[0], x_bh1[1], x_bh1[2], v_bh1[0], v_bh1[1], v_bh1[2], pot_act[0], -// m_bh2, x_bh2[0], x_bh2[1], x_bh2[2], v_bh2[0], v_bh2[1], v_bh2[2], pot_act[1]); -// fflush(stdout); - - - - for(i=2; i= t_stevol) - { - - /* Define the current mass of all star particles... */ - - for(i=0; i= t_stevol) */ - -#ifdef TIMING - get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst); - DT_STEVOL += (CPU_tmp_user - CPU_tmp_user0); -#endif - -#endif - - - /* STEVOL_SSE routine on all the nodes */ - -#ifdef STEVOL_SSE - - -#ifdef TIMING - get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); -#endif - - if(time_cur >= t_stevol) - { - - /* Define the current mass of all star particles... */ - - for(i=0; i= t_stevol) */ - -#ifdef TIMING - get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst); - DT_STEVOL += (CPU_tmp_user - CPU_tmp_user0); -#endif - - - - -/* Update ALL the new m,r,v "j" part. to the GRAPE/GPU memory */ - - - for(j=0; j= t_bh) - { - - if(myRank == rootRank) - { - -#ifdef BH_OUT - /* Write BH data... */ - write_bh_data(); -#endif - -#ifdef BH_OUT_NB - /* Write BH NB data... */ - write_bh_nb_data(); -#endif - -#ifdef BH_OUT_NB_EXT - /* Write BH NB data... */ - write_bh_nb_data_ext(); -#endif - - } /* if(myRank == rootRank) */ - - t_bh += dt_bh; - } /* if(time_cur >= t_bh) */ - - - - - - if(time_cur >= t_contr) - { - - if(myRank == rootRank) - { - -#ifdef STARDESTR - -/* - out = fopen("phi-GRAPE.ext","w"); - fprintf(out,"%.8E \t %.8E \n", m_bh, b_bh); - fclose(out); -*/ -/* - out = fopen("mass-bh.dat","a"); - fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); - fclose(out); -*/ - -#ifdef ADD_BH2 - - out = fopen("mass-bh.dat","a"); - fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", - time_cur, m_bh1, num_bh1, m_bh2, num_bh2); - fclose(out); - - out = fopen("mass-bh.con","w"); - fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", - time_cur, m_bh1, num_bh1, m_bh2, num_bh2); - fclose(out); - -#else - -#ifdef ADD_BH1 - - out = fopen("mass-bh.dat","a"); - fprintf(out,"%.8E \t %.8E %06d \n", - time_cur, m_bh1, num_bh1); - fclose(out); - - out = fopen("mass-bh.con","w"); - fprintf(out,"%.8E \t %.8E %06d \n", - time_cur, m_bh1, num_bh1); - fclose(out); - -#endif // ADD_BH1 - -#endif // ADD_BH2 - -#endif // STARDESTR - -#ifdef STARDESTR_EXT - - out = fopen("phi-GRAPE.ext","w"); - fprintf(out,"%.8E \t %.8E \n", m_bh, b_bh); - fclose(out); - - out = fopen("mass-bh.dat","a"); - fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); - fclose(out); - - out = fopen("mass-bh.con","w"); - fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); - fclose(out); - -#endif // STARDESTR_EXT - - - -#ifdef ACT_DEF_LL -#ifdef DEBUG111 - if(myRank == rootRank) check_linklist((int)Timesteps); -#endif -#endif - - - energy_contr(); - -#ifdef CMCORR111 - for(i=0; i R_LIMITS_ALL) && (m[i] != 0.0) ) - if( (tmp_r > R_LIMITS_ALL) && (tmp_rv > 0.0) ) - { - -// tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) ); -// E_corr += 0.5*m[i]*SQR(tmp_v); - -// pot[i] = 0.0; - -#ifdef EXTPOT -// pot_ext[i] = 0.0; -#endif - - for(k=0;k<3;k++) - { -// x[i][k] = R_LIMITS_ALL + 1.0*my_rand2(); - v[i][k] *= 1.0E-01; -// a[i][k] = 1.0E-06*my_rand2(); -// adot[i][k] = 1.0E-06*my_rand2(); - } /* k */ - -// t[i] = min_t; -// dt[i] = 0.125; - - } /* if */ - - } /* i */ - - - - for(j=0; j= t_contr) */ - - - - - - - - - if(time_cur >= t_disk) - { - - - if(myRank == rootRank) - { - - diskstep++; - - write_snap_data(); -// write_cont_data(); - -#ifdef DEBUG_extra - write_snap_extra(); -#endif - -#ifdef SPH_GAS - write_data_SPH(); -#endif - -#ifdef GMC - write_snap_GMC(); -#endif - - } /* if(myRank == rootRank) */ - -#ifdef ETICS_DUMP - sprintf(out_fname, "coeffs.%06d.%02d.dat", diskstep, myRank); - grapite_dump(out_fname, 2); -#endif - - - - -#ifdef LIMITS_NEW - - for(i=0; i 900.0 ) - { - - m[i] = 0.0; - - pot[i] = 0.0; - -#ifdef EXTPOT - pot_ext[i] = 0.0; -#endif - - for(k=0;k<3;k++) - { - x[i][k] = 1000.0 + 1.0*my_rand2(); - v[i][k] = 1.0E-06*my_rand2(); - a[i][k] = 1.0E-06*my_rand2(); - adot[i][k] = 1.0E-06*my_rand2(); - } /* k */ - - t[i] = 2.0*t_end; - - dt[i] = 0.125; - - } /* if */ - - } /* i */ - - - - for(j=0; j= t_disk) */ - - - -#ifdef CPU_TIMELIMIT - if(myRank == rootRank) - { - get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst); - tmp_cpu = CPU_time_real-CPU_time_real0; - } /* if(myRank == rootRank) */ -#endif - - -#ifdef ACT_DEF_LL - //printf("At tsteps= %04d . t+dt for BH = % 08E\n", int(Timesteps), t[N-1]+dt[N-1]); - ModifyLinkList(); - -#ifdef DEBUG111 - if( (((int)Timesteps)%10000 == 0) && (myRank == rootRank) ) check_linklist((int)Timesteps); -#endif -#endif - - } /* while(time_cur < t_end) */ - - - - /* close the local GRAPE's */ - -#ifdef SAPPORO - g6_close_(&clusterid); -#else - g6_close(clusterid); -#endif - - - /* Wait to all processors to finish his works... */ - MPI_Barrier(MPI_COMM_WORLD); - - MPI_Reduce(&g6_calls, &g6_calls_sum, 1, MPI_DOUBLE, MPI_SUM, rootRank, MPI_COMM_WORLD); - - /* Wait to all processors to finish his works... */ - MPI_Barrier(MPI_COMM_WORLD); - - - if(myRank == rootRank) - { - - /* Write some output for the timestep annalize... */ - - printf("\n"); - printf("Timesteps = %.0f Total sum of integrated part. = %.0f g6_calls on all nodes = %.0f \n", Timesteps, n_act_sum, g6_calls); - printf("\n"); - printf("Real Speed = %.3f GFlops \n", 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09); - fflush(stdout); - - -#ifdef DEBUG111 - tmp_file = fopen("n_act_distr.dat","w"); - - fprintf(tmp_file,"\n"); - fprintf(tmp_file,"Total Timesteps = %.0f Total sum of integrated part. = %.0f g6_calls on all nodes = %.0f \n", Timesteps, n_act_sum, g6_calls); - fprintf(tmp_file,"\n"); - fprintf(tmp_file,"Real Speed = %.3f GFlops \n", 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09); - - for(i=1;i Plummer + : Andreas+Fazeel... + : + : Stellar evolution added + : Stellar lifetimes: Raiteri, Villata & Navarro (1996) + : IMS mass loss: van den Hoeg & Groenewegen (1997) + : + : STARDESTR_EXT: Tidal disruption of stars by external BH... + : Chingis, Denis & Maxim... + : + : STARDESTR: Tidal disruption of stars by BH... + : Jose, Li Shuo & Shiyan Zhong + : + : STARDISK: Drag force... + : Chingis, Denis & Maxim... + : + : STARDISK: variable hz = HZ*(R/R_crit) up to R_crit... + : Taras, Andreas... + : + : Live BH (1 || 2) + ACC + EJECT... + : Li Shuo & Shiyan Zhong + : + : dt_min for BH (1 || 2)... + : + : added the PN calculus for the BBH + : PN0, PN1, PN2, PN2.5 (coded on the base of + : Gabor Kupi original routine) + : + : added the "name" array... + : + : added the GMC's calculus (GMC on CPU; GMC2 on GPU) + : for Alexey SC runs... and also for Fazeel Zurich runs... + : + : CPU_TIMELIMIT added for the Julich MW cluster runs... + : +Coded by : Peter Berczik +Version number : 19.04 +Last redaction : 2019.04.16 12:55 +*****************************************************************************/ +//#define CPU_RUNLIMIT 28000 // 24h = 24*60*60 = 86400 +#define CPU_RUNLIMIT 255000 // 3 days = 3*24*60*60 = 259200 +//#define CPU_RUNLIMIT 100 // 9.5 min +//#define CPU_RUNLIMIT 100 // 9.5 min + +//#define GMC // add the GMC's to the run +//#define GMC2 // add the GMC's to the run + +//#define NORM // Physical normalization +//#define STEVOL // Stellar evolution Do not tuch! Defined inside the Makefiles !!! +//#define STEVOL_SSE // Stellar evolution with SSE Do not tuch! Defined inside the Makefiles !!! + +//#define ADD_BH1 // add the Single BH + +#define ADD_BH2 // add the Binary BH's +#define ADD_N_BH // eps_BH = 0.0, but added only the Newtonian forces +// #define ADD_PN_BH // extra - added also the Post-Newton forces +// #define ADD_SPIN_BH // extra - added the SPIN for the BH's - DEFAULT !!! + +// #define BH_OUT // extra output for BH's (live) +// #define BH_OUT_NB // extra output for the BH's neighbours (live) + +// #define BBH_INF // BBH influence sphere... +// #define R_INF 10.0 // Factor for the influence sphere... if ( R < R_INF * DR_BBH ) +// #define R_INF2 (R_INF*R_INF) + +// #define DT_MIN_WARNING // dt < dt_min warning !!! + +//#define BH_OUT_NB_EXT // extra output for the BH's neighbours (external) + +//#define SAPPORO // Gaburov. Do not tuch! Defined inside the Makefiles !!! +//#define YEBISU // Nitadori. Do not tuch! Defined inside the Makefiles !!! +//#define GUINNESS // Nakasato. Do not tuch! Defined inside the Makefiles !!! + +//#define STARDESTR // disruption of stars by BH tidal forces +//#define R_TIDAL 1.0E-03 // Tidal radius of the BH's +//#define RMAX_OUT // extra data for def. r_max... + +//#define STARDESTR_EXT // disruption of stars by external BH tidal forces +//#define R_0 0.22 // Outer cut of the disk +//#define R_ACCR 0.0050 // Tidal Accr. rad of the BH's in units of R_0 +//#define R_TIDAL (R_0*R_ACCR) // Tidal radius of the BH's in NB units + +//#define STARDISK // analytic gas disk around BH + drag +//#define HZ (0.001*R_0) // Disk thickness... in NB units +//#define R_CRIT 0.0257314 // Critical radius of the disk (vert SG = vert BH force) +//#define R_CRIT 0.0 // i.e. hz=HZ=const... + +//#define SPH_GAS // SPH gas disk around BH + drag; experimental stuff !!! + +//#define EXTPOT // external potential (BH? or galactic?) +//#define EXTPOT_BH // BH - usually NB units + +//#define EXTPOT_GAL // Galactic B+D+H PK - usually physical units + +//#define EXTPOT_GAL_DEH // Dehnen Galactic - usually physical units +//#define EXTPOT_GAL_LOG // Log Galactic - usually physical units + +//#define EXTPOT_SC // SC extra POT for Bek test runs... +//#define DATAFILE eff0.05.tsf10.00.tab +//#define M_SC_DIM 100001 + +//#define CMCORR // CM correction in the zero step and in every dt_contr + +//#define DEBUG // debug information to files & screen +// #define DEBUG_extra // extra debug information to files & screen +//#define LAGR_RAD // write out lagr-rad.dat + +//#define LIMITS // for "0" mass particles !!! +//#define R_LIMITS 1.0E+03 // for "0" mass particles !!! +//#define COORDS 1010.0 // for "0" mass particles !!! + +//#define LIMITS_NEW // for "0" mass particles !!! + +//#define LIMITS_ALL_BEG // for ALL particles at the beginning... +//#define LIMITS_ALL // for ALL particles +//#define R_LIMITS_ALL 1.0E+03 // for ALL particles + +#ifdef ETICS +#include "grapite.h" +// why do we need CEP as a compilaion flag... just have it always on when ETICS is on. IF there is no CEP, there should be a graceful skipping of those operations. +//#define ETICS_CEP +#ifndef ETICS_DTSCF +#error "ETICS_DTSCF must be defined" +#endif +#endif + +#define TIMING + +#define ETA_S_CORR 4.0 +#define ETA_BH_CORR 4.0 + +#define DTMAXPOWER -3.0 +#define DTMINPOWER -36.0 + +/* +-3.0 0.125 +-4.0 0.0625 +-5.0 0.03125 +-7.0 ~1e-2 +-10.0 ~1e-3 +............. +-20.0 ~1e-6 +-23.0 ~1e-7 +-26.0 ~1e-8 +-30.0 ~1e-9 +*/ + +//#define AMD + +// #define ACT_DEF_LL + +#if defined(ACT_DEF_LL) && defined(ACT_DEF_GRAPITE) +#error "Contradicting preprocessor flags!" +#endif + +/****************************************************************************/ +#include +#include +#include +#include +#include +#include +#include + +/* +double aaa; +double aaapars[5]; + +extern void qwerty_(double *aaa, double *aaapars); +*/ + +#ifdef SAPPORO +#include +#include +#define GPU_PIPE 256 +#else +# ifdef YEBISU +# define G6_NPIPE 2048 +# else +# define G6_NPIPE 48 +# endif +#include "grape6.h" +#endif + +#ifdef GUINNESS +#define GPU_PIPE 128 +#endif + +#include + +/* Some "good" functions and constants... */ +#define SIG(x) ( ((x)<0) ? (-1):(1) ) +#define ABS(x) ( ((x)<0) ? (-x):(x) ) +#define MAX(a,b) ( ((a)>(b)) ? (a):(b) ) +#define MIN(a,b) ( ((a)<(b)) ? (a):(b) ) +#define SQR(x) ( (x)*(x) ) +#define POW3(x) ( (x)*SQR(x) ) + +#define Pi 3.14159265358979323846 +#define TWOPi 6.283185307179 +#define sqrt_TWOPi 2.506628274631 + +#ifdef NORM +//http://pdg.lbl.gov/2015/reviews/rpp2015-rev-astrophysical-constants.pdf + +#define G 6.67388E-11 // (m/s^2) * (m^2/kg) +#define Msol 1.988489E+30 // kg +#define Rsol 6.957E+08 // m +#define AU 149597870700.0 // m +#define pc 3.08567758149E+16 // m +#define Year 31556925.2 // s +#define c_feny 299792458.0 // m/s + +#define kpc (1.0E+03*pc) // m +#define km 1.0E+03 // km -> m +#define cm3 1.0E-06 // cm^3 -> m^3 +#define Myr (1.0E+06*Year) // s +#define Gyr (1.0E+09*Year) // s +#define R_gas 8.31447215 // J/(K*mol) +#define k_gas 1.380650424E-23 // J/K +#define N_A 6.022141510E+23 // 1/mol +#define mu 1.6605388628E-27 // kg +#define mp 1.67262163783E-27 // kg +#define me 9.1093821545E-31 // kg + +#define pc2 (pc*pc) +#define pc3 (pc*pc*pc) +#define kpc2 (kpc*kpc) +#define kpc3 (kpc*kpc*kpc) +#endif + +/* + 1KB = 1024 + 2KB = 2048 + 4KB = 4096 + 8KB = 8192 +16KB = 16384 +32KB = 32768 +64KB = 65536 +128KB = 131072 +256KB = 262144 +512KB = 524288 +1024KB = 1048576 -> 1MB +*/ + +#define KB 1024 +#define MB (KB*KB) + +#define N_MAX (6*MB) +#define N_MAX_loc (2*MB) + +//#define N_MAX (1200000) +//#define N_MAX_loc (1200000) + +//#define P_MAX 32 + +//#ifdef DEBUG + +/* extra code & data for DEBUG... */ + +#include "debug.h" + +//#ifdef DEBUG_extra +int ind_sort[N_MAX]; +double var_sort[N_MAX]; +//#endif + +//#endif + +#ifdef LAGR_RAD +int lagr_rad_N = 22; +double mass_frac[] = { 0.0001, 0.0003, 0.0005, 0.001, 0.003, 0.005, 0.01, 0.03, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 0.9999, 1.0 }; +double lagr_rad[22]; +double m_tot, m_cum, n_cum; +#endif + +int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank, + i, j, k, ni, nj, diskstep=0, power, jjj, iii, + skip_con=0, tmp_i; + +double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, + dt_bh, t_bh=0.0, dt_bh_tmp, + t_end, time_cur, dt_min, dt_max, min_t, min_t_loc, dt_new, min_dt, + eta_s, eta, eta_bh, + E_pot, E_pot_ext, E_kin, E_tot, E_tot_0, DE_tot, + E_tot_corr, E_tot_corr_0, DE_tot_corr, + E_tot_corr_sd, E_tot_corr_sd_0, DE_tot_corr_sd, + E_corr = 0.0, E_sd = 0.0, t_diss_on = 0.125, + e_kin_corr = 0.0, e_pot_corr = 0.0, e_pot_BH_corr = 0.0, + e_summ_corr = 0.0, + mcm, rcm_mod, vcm_mod, + rcm_sum=0.0, vcm_sum=0.0, + eps=0.0, eps2, + xcm[3], vcm[3], mom[3], + xdc[3], vdc[2], + over2=(1.0/2.0), over3=(1.0/3.0), over6=(1.0/6.0), + a2_mod, adot2_mod, + dt_tmp, dt2half, dt3over6, dt4over24, dt5over120, + dtinv, dt2inv, dt3inv, + a0mia1, ad04plad12, ad0plad1, + a2[3], a3[3], a2dot1[3], + a1abs, adot1abs, a2dot1abs, a3dot1abs, + Timesteps=0.0, n_act_sum=0.0, n_act_distr[N_MAX], g6_calls=0.0, g6_calls_sum=0.0, + tmp, tmp_r, tmp_v, tmp_rv, tmp_cpu, + tmp_pot, tmp_a, tmp_adot, + tmp_a_bh, tmp_adot_bh, + tmp_a_bh_pn0, tmp_a_bh_pn1, tmp_a_bh_pn2, tmp_a_bh_pn2_5, tmp_a_bh_pn3, tmp_a_bh_pn3_5, tmp_a_bh_spin, + tmp_adot_bh_pn0, tmp_adot_bh_pn1, tmp_adot_bh_pn2, tmp_adot_bh_pn2_5, tmp_adot_bh_pn3, tmp_adot_bh_pn3_5, tmp_adot_bh_spin; + +double R[3], V[3], L[3], mju, dr2, dr, dv2, dv, dl2, dl, pot_eff; + +char processor_name[MPI_MAX_PROCESSOR_NAME], + inp_fname[30], inp_fname_gmc[30], + out_fname[30], dbg_fname[30]; + +/* global variables */ + +int N, N_star, N_gmc, N_bh, + ind[N_MAX], name[N_MAX]; +double m[N_MAX], x[N_MAX][3], v[N_MAX][3], + pot[N_MAX], a[N_MAX][3], adot[N_MAX][3], + t[N_MAX], dt[N_MAX]; + +/* local variables */ + +int n_loc, ind_loc[N_MAX_loc]; +double m_loc[N_MAX_loc], x_loc[N_MAX_loc][3], v_loc[N_MAX_loc][3], + pot_loc[N_MAX_loc], a_loc[N_MAX_loc][3], adot_loc[N_MAX_loc][3], + t_loc[N_MAX_loc], dt_loc[N_MAX_loc]; + +/* data for active particles */ + +int n_act, ind_act[N_MAX]; +double m_act[N_MAX], + x_act[N_MAX][3], v_act[N_MAX][3], + pot_act[N_MAX], a_act[N_MAX][3], adot_act[N_MAX][3], + t_act[N_MAX], dt_act[N_MAX], + x_act_new[N_MAX][3], v_act_new[N_MAX][3], + pot_act_new[N_MAX], a_act_new[N_MAX][3], adot_act_new[N_MAX][3], + pot_act_tmp[N_MAX], a_act_tmp[N_MAX][3], adot_act_tmp[N_MAX][3], + pot_act_tmp_loc[N_MAX], a_act_tmp_loc[N_MAX][3], adot_act_tmp_loc[N_MAX][3]; + +FILE *inp, *out, *tmp_file, *dbg; + +double CPU_time_real0, CPU_time_user0, CPU_time_syst0; +double CPU_time_real, CPU_time_user, CPU_time_syst; + + +#ifdef TIMING + +double CPU_tmp_real0, CPU_tmp_user0, CPU_tmp_syst0; +double CPU_tmp_real, CPU_tmp_user, CPU_tmp_syst; + +double DT_TOT, + DT_ACT_DEF1, DT_ACT_DEF2, DT_ACT_DEF3, DT_ACT_PRED, + DT_ACT_GRAV, DT_EXT_GRAV, + DT_GMC_GRAV, DT_GMC_GMC_GRAV, DT_EXT_GMC_GRAV, + DT_ACT_CORR, DT_ACT_LOAD, + DT_STEVOL, DT_STARDISK, DT_STARDESTR; + +double DT_ACT_REDUCE; + +#endif + +/* some local settings for G6a board's */ + +int clusterid, ii, nn, numGPU; + +//#ifdef SAPPORO +#if defined(SAPPORO) || defined(GUINNESS) + +int npipe=GPU_PIPE, index_i[GPU_PIPE]; +double h2_i[GPU_PIPE], x_i[GPU_PIPE][3], v_i[GPU_PIPE][3], + p_i[GPU_PIPE], a_i[GPU_PIPE][3], jerk_i[GPU_PIPE][3]; +double new_tunit=51.0, new_xunit=51.0; + +#elif defined YEBISU + +int npipe=G6_NPIPE, index_i[G6_NPIPE]; +double h2_i[G6_NPIPE], x_i[G6_NPIPE][3], v_i[G6_NPIPE][3], + p_i[G6_NPIPE], a_i[G6_NPIPE][3], jerk_i[G6_NPIPE][3]; +int new_tunit=51, new_xunit=51; +#else + +int npipe=48, index_i[48]; +double h2_i[48], x_i[48][3], v_i[48][3], + p_i[48], a_i[48][3], jerk_i[48][3]; +int new_tunit=51, new_xunit=51; + +#endif + +int aflag=1, jflag=1, pflag=1; + +double ti=0.0, a2by18[3], a1by6[3], aby2[3]; + +/* normalization... */ + +#ifdef NORM +double m_norm, r_norm, v_norm, t_norm; +#endif + +double eps_BH=0.0; + +/* external potential... */ + +#ifdef EXTPOT + +#ifdef EXTPOT_GAL +double m_bulge, a_bulge, b_bulge, + m_disk, a_disk, b_disk, + m_halo, a_halo, b_halo, +// x_ext, y_ext, z_ext, +// vx_ext, vy_ext, vz_ext, +// x_ij, y_ij, z_ij, +// vx_ij, vy_ij, vz_ij, rv_ij, + x2_ij, y2_ij, z2_ij, + r_tmp, r2_tmp, z_tmp, z2_tmp; +#endif + +#ifdef EXTPOT_GAL_DEH +double m_ext, r_ext, g_ext, + tmp_r2, tmp_r3, dum, dum2, dum3, dum_g, tmp_g; +#endif + +#ifdef EXTPOT_GAL_LOG +double v_halo, r_halo, + v2_halo, r2_halo, r2_r2_halo, + x2_ij, y2_ij, z2_ij; +#endif + +#ifdef EXTPOT_BH +double r2, rv_ij, + m_bh, b_bh, eps_bh; +#endif + +#ifdef EXTPOT_SC +int M_SC=M_SC_DIM; +double r_sc[M_SC_DIM], m_sc[M_SC_DIM], p_sc[M_SC_DIM]; +double M_R, pot_out_R; +#endif + + +double pot_ext[N_MAX], pot_act_ext[N_MAX]; // for all EXTPOT + +#endif + + +int i_bh, i_bh1, i_bh2, + num_bh = 0, num_bh1 = 0, num_bh2 = 0; + +double m_bh, m_bh1, m_bh2, b_bh, + r, r2, + x_ij, y_ij, z_ij, + vx_ij, vy_ij, vz_ij, rv_ij; + +#ifdef STEVOL +int num_dead, event[N_MAX]; +double dt_stevol, t_stevol; +double m0[N_MAX], ZZZ[N_MAX], t0[N_MAX]; +#endif + + +#ifdef STEVOL_SSE + +int num_dead, event[N_MAX], event_old[N_MAX]; +double dt_stevol, t_stevol; +double m0[N_MAX], ZZZ[N_MAX], t0[N_MAX]; +double SSEMass[N_MAX], SSESpin[N_MAX], SSERad[N_MAX], + SSELum[N_MAX], SSETemp[N_MAX]; +double SSEMV[N_MAX], SSEBV[N_MAX]; + +double Rgal = 10000.0; //Distance of star from sun for artificial CMD with observational errors [pc] + +int van_kick=0; + +extern void zcnsts_(double *z, double *zpars); +extern void evolv1_(int *kw, double *mass, double *mt, double *r, double *lum, double *mc, double *rc, double *menv, double *renv, double *ospin, double *epoch, double *tms, double *tphys, double *tphysf, double *dtp, double *z, double *zpars, double *vkick, double *vs); +extern int cmd(int kstar, double rstar, double lstar, double Rgal, double *abvmag, double *vmag, double *BV, double *Teff, double *dvmag, double *dBV); + +struct{ + double neta; + double bwind; + double hewind; + double mxns; +} value1_; + +struct{ + double pts1; + double pts2; + double pts3; +} points_; + +struct{ + double sigma; + int bhflag; +} value4_; + +struct{ + int idum; +} value3_; + +struct{ + int ceflag; + int tflag; + int ifflag; + int nsflag; + int wdflag; +} flags_; + +struct{ + double beta; + double xi; + double acc2; + double epsnov; + double eddfac; + double gamma; +} value5_; + +struct{ + double alpha1; + double lambda; +} value2_; + +#include "cmd.c" + +#endif + + + +#ifdef BBH_INF +int inf_event[N_MAX]; +double x_bbhc[3], v_bbhc[3], DR2, tmp_r2; +double DV2, EB, SEMI_a, SEMI_a2; +#endif + + + +/* STARDISK Drag force... */ + +#ifdef STARDISK + +double a_drag[N_MAX][3], adot_drag[N_MAX][3]; +//double a_drag[N_MAX][3], adot_drag[N_MAX][3]; + +double z_new_drag, r_new_drag2, hz; + +#ifdef SPH_GAS + +#define N_GAS_MAX (128*KB) +#define NB 50 + +int N_GAS, ind_gas[N_GAS_MAX], sosed_gas[N_GAS_MAX][NB]; +double m_gas[N_GAS_MAX], x_gas[N_GAS_MAX][3], v_gas[N_GAS_MAX][3], h_gas[N_GAS_MAX], DEN_gas[N_GAS_MAX]; + +//int sosed[NB]; +double d[N_MAX]; // podrazumevajem chto: N_MAX > N_GAS_MAX vsegda ! +double x_star, y_star, z_star, DEN_star; +double h_min = 1.0E-04, h_max = 1.0; + +#include "def_H.c" +#include "def_DEN.c" +#include "def_DEN_STAR.c" + +#endif // SPH_GAS + +#include "drag_force.c" + +// calculate the SG for the set of active particles which is deepley INSIDE the EXT BH infl. rad. +// EXPERIMENTAL feature !!! + +//int SG_CALC = 0; +//double summ_tmp_r = 0.0, aver_tmp_r = 0.0, R_INT_CUT = 1.0E-03; + +#endif // STARDISK + + + +/* GMC data... */ + +#ifdef GMC + +double dt_gmc, E_pot_gmc, E_pot_gmc_gmc, E_pot_ext_gmc, E_kin_gmc; + +#define N_gmc_MAX (500) + +int N_gmc, ind_gmc[N_gmc_MAX], name_gmc[N_gmc_MAX]; + +double m_gmc[N_gmc_MAX], pot_gmc_gmc[N_gmc_MAX], pot_ext_gmc[N_gmc_MAX], + x_gmc[N_gmc_MAX][3], v_gmc[N_gmc_MAX][3], a_gmc[N_gmc_MAX][3], + eps_gmc[N_gmc_MAX]; + +double pot_gmc[N_MAX]; + +#endif // GMC + + +/****************************************************************************/ + +/****************************************************************************/ +#ifdef ADD_N_BH + +double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3]; + +double pot_bh1, a_bh1[3], adot_bh1[3], + pot_bh2, a_bh2[3], adot_bh2[3]; + +//double eps_BH = 0.0; + +#include "n_bh.c" + +/* +int calc_force_n_BH(double m1, double xx1[], double vv1[], + double m2, double xx2[], double vv2[], + double eps_BH, + double pot_n1, double a_n1[], double adot_n1[], + double pot_n2, double a_n2[], double adot_n2[]) +*/ + +/* + INPUT + +m1 - mass of the 1 BH +xx1[0,1,2] - coordinate of the 1 BH +vv1[0,1,2] - velocity of the 1 BH + +m2 - mass of the 2 BH +xx2[0,1,2] - coordinate of the 2 BH +vv2[0,1,2] - velocity of the 2 BH + +eps_BH - force softening, can be even exactly 0.0 ! + + OUTPUT + +pot_n1 for the 1 BH +a_n1 [0,1,2] for the 1 BH +adot_n1 [0,1,2] for the 1 BH + +pot_n2 for the 2 BH +a_n2 [0,1,2] for the 2 BH +adot_n2 [0,1,2] for the 2 BH + +return - 0 if everything OK +*/ + +#endif +/****************************************************************************/ + +/****************************************************************************/ +#ifdef ADD_PN_BH + +double C_NB = 477.12; + +int usedOrNot[7] = {1, 1, 1, 1, 0, 0, 0}; + +double a_pn1[7][3], adot_pn1[7][3], + a_pn2[7][3], adot_pn2[7][3]; + +double s_bh1[3] = {0.0, 0.0, 1.0}; +double s_bh2[3] = {0.0, 0.0, 1.0}; + +#ifdef ADD_SPIN_BH // eto rabotajet vsegda !!! +#include "pn_bh_spin.c" +#else +#include "pn_bh.c" // eto staraja versija, bolshe ne rabotajet !!! +#endif + +/* +int calc_force_pn_BH(double m1, double xx1[], double vv1[], double ss1[], + double m2, double xx2[], double vv2[], double ss2[], + double CCC_NB, double dt_bh, + int usedOrNot[], + double a_pn1[][3], double adot_pn1[][3], + double a_pn2[][3], double adot_pn2[][3]) +*/ + +/* + INPUT + +m1 - mass of the 1 BH +xx1[0,1,2] - coordinate of the 1 BH +vv1[0,1,2] - velocity of the 1 BH +spin1[0,1,2] - normalized spin of the 1 BH + +m2 - mass of the 2 BH +xx2[0,1,2] - coordinate of the 2 BH +vv2[0,1,2] - velocity of the 2 BH +spin2[0,1,2] - normalized spin of the 2 BH + +CCC_NB - Speed of light "c" in internal units +dt_BH - timestep of the BH's, needed for the SPIN integration + +usedOrNot[PN0, PN1, PN2, PN2.5, PN3, PN3.5, SPIN] - different PN term usage: PN1, PN2, PN2.5, PN3, PN3.5, SPIN + 0 1 2 3 4 5 6 + + OUTPUT + +a_pn1 [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 1 BH +adot_pn1[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 1 BH + +a_pn2 [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 2 BH +adot_pn2[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 2 BH + +return - 0 if everything OK + - 505 if BH's separation < 4 x (RSwarch1 + RSwarch2) +*/ + +#endif + +/****************************************************************************/ + +/****************************************************************************/ +/* RAND_MAX = 2147483647 */ +/* my_rand : 0.0 - 1.0 */ +/* my_rand2 : -1.0 - 1.0 */ + +double my_rand(void) +{ +return( (double)(rand()/(double)RAND_MAX) ); +} + +double my_rand2(void) +{ +return (double)(2.0)*((rand() - RAND_MAX/2)/(double)RAND_MAX); +} +/****************************************************************************/ + +#ifdef STARDESTR +#include "star_destr.c" +#endif + +#ifdef STARDESTR_EXT +#include "star_destr_ext.c" +#endif + +#ifdef ACT_DEF_LL +#include "act_def_linklist.c" +#endif + +#ifdef ETICS +double t_exp, dt_exp=ETICS_DTSCF; // t_exp is just the initial value +#ifdef ETICS_CEP +int grapite_cep_index; +#endif +#endif + +/****************************************************************************/ +void get_CPU_time(double *time_real, double *time_user, double *time_syst) +{ +struct rusage xxx; +double sec_u, microsec_u, sec_s, microsec_s; +struct timeval tv; + + +getrusage(RUSAGE_SELF,&xxx); + +sec_u = xxx.ru_utime.tv_sec; +sec_s = xxx.ru_stime.tv_sec; + +microsec_u = xxx.ru_utime.tv_usec; +microsec_s = xxx.ru_stime.tv_usec; + +*time_user = sec_u + microsec_u * 1.0E-06; +*time_syst = sec_s + microsec_s * 1.0E-06; + +// *time_real = time(NULL); + +gettimeofday(&tv, NULL); +*time_real = tv.tv_sec + 1.0E-06 * tv.tv_usec; + +*time_user = *time_real; + +} +/****************************************************************************/ + +#ifdef STEVOL +/****************************************************************************/ +double t_dead(double m, double Z) +/* +m - mass of the star (in Msol) +Z - metallicity (absolut values - i.e. Zsol = 0.0122) +t - return value of star lifetime (in Year) + +Raiteri, Villata & Navarro, 1996, A&A, 315, 105 +*/ +{ +double a0, a1, a2, lZ, lm, tmp=0.0; + +m *= (m_norm/Msol); + +lZ = log10(Z); +lm = log10(m); + +a0 = 10.130 + 0.07547*lZ - 0.008084*lZ*lZ; +a1 = -4.424 - 0.79390*lZ - 0.118700*lZ*lZ; +a2 = 1.262 + 0.33850*lZ + 0.054170*lZ*lZ; + +tmp = a0 + a1*lm + a2*lm*lm; + +tmp = pow(10.0,tmp); + +tmp *= (Year/t_norm); + +return(tmp); +} +/****************************************************************************/ + +/****************************************************************************/ +double m_loss_old(double m, double Z) +/* +m - initial mass of the star (in Msol) +Z - metallicity (e.g. Zsol = 0.02) +m_ret - returned maass to ISM from the star (in Msol) + +approx. from data of van den Hoek & Groenewegen, 1997, A&AS, 123, 305 +*/ +{ +double tmp=0.0; + +m *= (m_norm/Msol); + +// tmp = (-0.46168 + 0.912274 * m); // first approx. + +tmp = -0.420542*pow(Z, -0.0176601) + (0.901495 + 0.629441*Z) * m; + +tmp *= (Msol/m_norm); + +return(tmp); +} +/****************************************************************************/ + +/****************************************************************************/ +double m_curr_old(double m0, double Z, double t0, double t) +{ +double td, ml, tmp=0.0; + +td = t_dead(m0, Z); +ml = m_loss_old(m0, Z); + +if ( (t-t0) < td ) + { + tmp = m0; // instant mass loss +// tmp = m0 - (t-t0)/(td-t0) * ml; // linear mass loss + } +else + { + tmp = m0 - ml; + + if (event[i] == 0) + { + num_dead++; + event[i] = 1; + } + + } /* if ( (t-t0) < td ) */ + +return(tmp); +} +/****************************************************************************/ + +/****************************************************************************/ +double m_loss(double m, double Z) +/* +m - initial mass of the star (in norm. units) +Z - metallicity (absolut value - i.e. Zsol = 0.0122) +m_loss - returned maass to ISM from the star (in norm. units) +*/ +{ +double tmp=0.0; + +m *= (m_norm/Msol); + +if ( m < 8.0) + { +// approx. from data of van den Hoek & Groenewegen, 1997, A&AS, 123, 305 + +// tmp = (-0.46168 + 0.912274 * m); // first approx. + tmp = -0.420542*pow(Z, -0.0176601) + (0.901495 + 0.629441*Z) * m; // second approx. + } +else + { +// approx. from data of Woosley & Weaver, 1995, ApJS, 110, 181 + +// tmp = 0.813956*pow(m, 1.04214); + tmp = 0.813956*pow(m, 1.04); + } + +tmp *= (Msol/m_norm); + +return(tmp); +} +/****************************************************************************/ + +/****************************************************************************/ +double m_curr(double m0, double Z, double t0, double t) +{ +double td, ml, tmp=0.0; + +tmp = m0; + +if ( event[i] == 0 ) + { + + td = t_dead(m0, Z); + + if ( (t-t0) > td ) + { + + ml = m_loss(m0, Z); + + tmp = m0 - ml; + + num_dead++; + +// wd http://en.wikipedia.org/wiki/Chandrasekhar_limit 1.38 + + if ( tmp * (m_norm/Msol) < 1.38 ) + { + event[i] = 1; + } + else + { + +// ns http://en.wikipedia.org/wiki/Neutron_stars 1.38 - ~2.0 ??? + + if ( tmp * (m_norm/Msol) > 1.38 ) event[i] = 2; + +// bh http://en.wikipedia.org/wiki/Tolman-Oppenheimer-Volkoff_limit ~1.5 - ~3.0 ??? + + if ( tmp * (m_norm/Msol) > 2.00 ) event[i] = 3; + + } + + } /* if ( (t-t0) > td ) */ + + } /* if ( event[i] == 0 ) */ + +return(tmp); +} +/****************************************************************************/ +#endif + + +#ifdef STEVOL_SSE +/****************************************************************************/ +double m_curr(int ipart, double m0, double Z, double t0, double t) +{ +//set up parameters and variables for SSE (Hurley, Pols & Tout 2002) +int kw=0; //stellar type +double mass; //initial mass +double mt=0.0; //actual mass +double r=0.0; //radius +double lum=0.0; //luminosity +double mc=0.0; //core mass +double rc=0.0; //core radius +double menv=0.0; //envelope mass +double renv=0.0; //envelope radius +double ospin=0.0; //spin +double epoch=0.0; //time spent in current evolutionary state +double tms=0.0; //main-sequence lifetime +double tphys; //initial age +double tphysf; //final age +double dtp=0.0; //data store value, if dtp>tphys no data will be stored +double z; //metallicity +double zpars[20]; //metallicity parameters +double vkick=0.0; //kick velocity for compact remnants +double vs[3]; //kick velocity componets + +// M_V_[mag] V_[mag] B-V_[mag] T_eff_[K] dV_[mag] d(B-V) +double abvmag, vmag, BV, Teff, dvmag, dBV; + +mass = m0*(m_norm/Msol); // m0[i] initial mass of the stars Msol +mt = mass; // current mass for output in Msol +tphys = t0*(t_norm/Myr); // t0[i] time of the star formation Myr +tphysf = t*(t_norm/Myr); // t[i] current time for the star Myr +dtp = 0.0; // output parameter, just set to 0.0 !!! +z = Z; // ZZZ[i] of the star + +/* +if (ipart > 9995) + { +printf("Initial parameters: \n"); +printf("M = %g [Mo] \n", mass); +printf("Z = %g \n", z); +printf("kw = %d \n", kw); +printf("spin = %g \n", ospin); +printf("epoch = %g [Myr] \n", epoch); +printf("t_beg = %g [Myr] \n", tphys); +printf("t_end = %g [Myr] \n", tphysf); +printf("R = %g [R_o] \n", r); +printf("L = %g [L_o] \n", lum); +printf("V_kick = %g [km/s] \n", vkick); + } +*/ + +for (k=0; k<20; k++) zpars[k] = 0; +zcnsts_(&z,zpars); //get metallicity parameters + +evolv1_(&kw, &mass, &mt, &r, &lum, &mc, &rc, &menv, &renv, &ospin, &epoch, &tms, &tphys, &tphysf, &dtp, &z, zpars, &vkick, vs); +cmd(kw, r, lum, Rgal, &abvmag, &vmag, &BV, &Teff, &dvmag, &dBV); + +/* +if (ipart > 9995) + { +printf("Final parameters: \n"); +printf("M = %g [Mo] \n", mt); +printf("Z = %g \n", z); +printf("kw = %d \n", kw); +printf("spin = %g \n", ospin); +printf("epoch = %g [Myr] \n", epoch); +printf("t_beg = %g [Myr] \n", tphys); +printf("t_end = %g [Myr] \n", tphysf); +printf("R = %g [R_o] \n", r); +printf("L = %g [L_o] \n", lum); +printf("V_kick = %g [km/s] \n", vkick); +printf("M_V_[mag]\tV_[mag]\t\tB-V_[mag]\tT_eff_[K]\tdV_[mag]\td(B-V)\n"); +printf("%.3E\t%.3E\t%.3E\t%.3E\t%.3E\t%.3E\n", abvmag, vmag, BV, Teff, dvmag, dBV); +// printf("%g\t%g\t%g\t%g\t%g\t%g\n", abvmag, vmag, BV, Teff, dvmag, dBV); + } +*/ + +event[ipart] = kw; // Evolution type (0 - 15). +SSEMass[ipart] = mt; // Mass in Msol +SSESpin[ipart] = ospin; // Spin in SI ??? [kg*m*m/s] +SSERad[ipart] = r; // Radius in Rsol +SSELum[ipart] = lum; // Luminosity in Lsol +SSETemp[ipart] = Teff; // Temperature + +SSEMV[ipart] = abvmag; // M_V absolute V magnitude +SSEBV[ipart] = BV; // B-V color index + +/* +c STELLAR TYPES - KW +c +c 0 - deeply or fully convective low mass MS star +c 1 - Main Sequence star +c 2 - Hertzsprung Gap +c 3 - First Giant Branch +c 4 - Core Helium Burning +c 5 - First Asymptotic Giant Branch +c 6 - Second Asymptotic Giant Branch +c 7 - Main Sequence Naked Helium star +c 8 - Hertzsprung Gap Naked Helium star +c 9 - Giant Branch Naked Helium star +c 10 - Helium White Dwarf +c 11 - Carbon/Oxygen White Dwarf +c 12 - Oxygen/Neon White Dwarf +c 13 - Neutron Star +c 14 - Black Hole +c 15 - Massless Supernova +*/ + +if ( van_kick == 1 ) // we have a kick ??? + { + +if ( (event_old[ipart] < 10) && ( (event[ipart] == 13) || (event[ipart] == 14) ) ) // NS or BH + { + + if (myRank == rootRank) + { +// printf("KICK: %06d %.6E [Myr] %02d (%02d) %.6E [Mo] %.6E [km/s] [% .3E, % .3E, % .3E] OLD: [% .3E, % .3E, % .3E] %.6E [Mo] \n", ipart, tphysf, kw, event_old[ipart], mt, vkick, vs[0], vs[1], vs[2], v[ipart][0]*(v_norm/km), v[ipart][1]*(v_norm/km), v[ipart][2]*(v_norm/km), mass); + printf("KICK: %06d %.6E %02d %02d %.6E %.6E % .3E % .3E % .3E % .3E % .3E % .3E %.6E \n", ipart, tphysf, kw, event_old[ipart], mt, vkick, vs[0], vs[1], vs[2], v[ipart][0]*(v_norm/km), v[ipart][1]*(v_norm/km), v[ipart][2]*(v_norm/km), mass); + fflush(stdout); + } /* if (myRank == rootRank) */ + + v[ipart][0] += vs[0]*(km/v_norm); + v[ipart][1] += vs[1]*(km/v_norm); + v[ipart][2] += vs[2]*(km/v_norm); + } + + } // we have a kick ??? + + +event_old[ipart] = event[ipart]; + + +tmp = mt*(Msol/m_norm); // return the current mass(t-t0) of the star in NB units + +return(tmp); +} +/****************************************************************************/ +#endif + + +/****************************************************************************/ +void read_data() +{ +inp = fopen(inp_fname,"r"); +fscanf(inp,"%d \n", &diskstep); + +//#ifdef STEVOL +// fscanf(inp,"%d \n", &N); +//#endif +// fscanf(inp,"%d %d %d %d \n", &N, &N_bh, &N_gmc, &N_star); + +fscanf(inp,"%d \n", &N); + +fscanf(inp,"%lE \n", &time_cur); + +for (i=0; i= (mass_frac[k]*m_tot) ) + { + lagr_rad[k] = var_sort[j]; + k++; + } +#endif + } + +fclose(out); + +// write out lagr-rad.dat + +#ifdef LAGR_RAD +out = fopen("lagr-rad.dat","a"); +fprintf(out,"%.6E ", time_cur); +for (k=0; k 100.0) jerk_i[ii][0] = 100.0; + if (ABS(jerk_i[ii][1]) > 100.0) jerk_i[ii][1] = 100.0; + if (ABS(jerk_i[ii][2]) > 100.0) jerk_i[ii][2] = 100.0; + } + +#ifdef SAPPORO + g6calc_firsthalf_(&clusterid, &n_loc, &nn, index_i, x_i, v_i , a_i, jerk_i, p_i, &eps2, h2_i); + g6calc_lasthalf_(&clusterid, &n_loc, &nn, index_i, x_i, v_i, &eps2, h2_i, a_i, jerk_i, p_i); +#else + g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i); + g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i); +#endif + + for (ii=0; ii r) ) + { + r1 = r_sc[ii]; r2 = r_sc[ii+1]; + m1 = m_sc[ii]; m2 = m_sc[ii+1]; + p1 = p_sc[ii]; p2 = m_sc[ii+1]; + break; + } + } + +m_tmp = m1 + (r-r1)*(m2-m1)/(r2-r1); +p_tmp = p1 + (r-r1)*(p2-p1)/(r2-r1); + +// if (myRank == rootRank) +// { +// printf("\t \t \t %.6E %.6E \n", r1, r2); +// printf("\t \t \t %.6E %.6E \n", m1, m2); +// printf("\t \t \t %.6E %.6E \n", p1, p2); +// printf("\t \t \t %06d %.6E %.6E %.6E \n", ii, r, m_tmp, p_tmp); +// fflush(stdout); +// } /* if (myRank == rootRank) */ + +*M_R = m_tmp; +*pot_out_R = p_tmp; + +//exit(-1); + +} +#endif +/****************************************************************************/ + +/****************************************************************************/ +void calc_ext_grav_zero() +{ + +#ifdef EXTPOT + +/* Define the external potential for all particles on all nodes */ + +ni = n_act; + +for (i=0; i 100.0) jerk_i[ii][0] = 100.0; + if (ABS(jerk_i[ii][1]) > 100.0) jerk_i[ii][1] = 100.0; + if (ABS(jerk_i[ii][2]) > 100.0) jerk_i[ii][2] = 100.0; + } + + g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i); + g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i); + + g6_calls++; +*/ + + for (ii=0; ii R_LIMITS_ALL) && (tmp_rv > 0.0) ) + { + + for (k=0;k<3;k++) + { + v[i][k] *= 1.0E-01; + } /* k */ + + } /* if */ + + } /* i */ + + #endif + + + #ifdef CMCORR + + /* possible CM correction of the initial datafile... */ + + if ( (diskstep == 0) && (time_cur == 0.0) ) cm_corr(); + + #endif + + + printf("\n"); + printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc); + printf("\n"); + #ifdef GMC + printf("N = %06d \t N_GMC = %06d \t eps = %.6E \n", N, N_gmc, eps); + #else + printf("N = %07d \t eps = %.6E \n", N, eps); + #endif + printf("t_beg = %.6E \t t_end = %.6E \n", time_cur, t_end); + printf("dt_disk = %.6E \t dt_contr = %.6E \n", dt_disk, dt_contr); + printf("dt_bh = %.6E \n", dt_bh); + printf("eta = %.6E \n", eta); + printf("\n"); + + #ifdef NORM + printf("Normalization: \n"); + printf("\n"); + printf("m_norm = %.4E [Msol] r_norm = %.4E [pc] \n", m_norm/Msol, r_norm/pc); + printf("v_morm = %.4E [km/s] t_morm = %.4E [Myr] \n", v_norm/km, t_norm/Myr); + printf("\n"); + #endif + + #ifdef EXTPOT + printf("External Potential: \n"); + printf("\n"); + + #ifdef EXTPOT_GAL + printf("m_bulge = %.4E [Msol] a_bulge = %.4E b_bulge = %.4E [kpc] \n", m_bulge*m_norm/Msol, a_bulge*r_norm/kpc, b_bulge*r_norm/kpc); + printf("m_disk = %.4E [Msol] a_disk = %.4E b_disk = %.4E [kpc] \n", m_disk*m_norm/Msol, a_disk*r_norm/kpc, b_disk*r_norm/kpc); + printf("m_halo = %.4E [Msol] a_halo = %.4E b_halo = %.4E [kpc] \n", m_halo*m_norm/Msol, a_halo*r_norm/kpc, b_halo*r_norm/kpc); + #endif + + #ifdef EXTPOT_BH + printf("m_bh = %.4E b_bh = %.4E \n", m_bh, b_bh); + #endif + + #ifdef EXTPOT_GAL_DEH + printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n", m_ext, r_ext, g_ext); + #endif + + #ifdef EXTPOT_GAL_LOG + printf("v_halo = %.6E r_halo = %.6E \n", v_halo, r_halo); + #endif + + #ifdef EXTPOT_SC + read_SC_mass(); + + printf("EXTPOT_SC: # of points %06d \t M_SC_TOT = %.6E \n", M_SC_DIM-1, m_sc[M_SC_DIM-1]); + #endif + + printf("\n"); + #endif + + + + fflush(stdout); + + + if ( (diskstep == 0) && (time_cur == 0.0) ) + { + + #ifdef SPH_GAS + write_data_SPH(); + #endif + + #ifdef GMC + write_snap_GMC(); + write_cont_GMC(); + #endif + + + out = fopen("contr.dat","w"); + fclose(out); + + #ifdef TIMING + out = fopen("timing.dat","w"); + fclose(out); + #endif + + #ifdef ADD_BH2 + + #ifdef BH_OUT + out = fopen("bh.dat","w"); + fclose(out); + #endif + + #ifdef BH_OUT_NB + out = fopen("bh_nb.dat","w"); + fclose(out); + #endif + + #else + + #ifdef ADD_BH1 + + #ifdef BH_OUT + out = fopen("bh.dat","w"); + fclose(out); + #endif + + #ifdef BH_OUT_NB + out = fopen("bh_nb.dat","w"); + fclose(out); + #endif + + #endif // ADD_BH1 + + #endif // ADD_BH2 + + + #ifdef STARDESTR + + #ifdef ADD_BH2 + + out = fopen("mass-bh.dat","w"); + m_bh1 = m[0]; + m_bh2 = m[1]; + num_bh1 = 0; + num_bh2 = 0; + fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", + time_cur, m_bh1, num_bh1, m_bh2, num_bh2); + fclose(out); + + out = fopen("mass-bh.con","w"); + m_bh1 = m[0]; + m_bh2 = m[1]; + num_bh1 = 0; + num_bh2 = 0; + fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", + time_cur, m_bh1, num_bh1, m_bh2, num_bh2); + fclose(out); + + #else + + #ifdef ADD_BH1 + + out = fopen("mass-bh.dat","w"); + m_bh1 = m[0]; + num_bh1 = 0; + fprintf(out,"%.8E \t %.8E %06d \n", + time_cur, m_bh1, num_bh1); + fclose(out); + + out = fopen("mass-bh.con","w"); + m_bh1 = m[0]; + num_bh1 = 0; + fprintf(out,"%.8E \t %.8E %06d \n", + time_cur, m_bh1, num_bh1); + fclose(out); + + #endif // ADD_BH1 + + #endif // ADD_BH2 + + #endif // STARDESTR + + + #ifdef STARDESTR_EXT + + num_bh = 0; + out = fopen("mass-bh.dat","w"); + fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); + fclose(out); + + num_bh = 0; + out = fopen("mass-bh.con","w"); + fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); + fclose(out); + + #ifdef BH_OUT_NB_EXT + out = fopen("bh_nb.dat","w"); + fclose(out); + #endif + + #endif // STARDESTR_EXT + + + #ifdef DEBUG_extra + #ifdef LAGR_RAD + out = fopen("lagr-rad.dat","w"); + fclose(out); + #endif + #endif + + + #ifdef BBH_INF + out = fopen("bbh.inf","w"); + fclose(out); + #endif + + + } + else // if (diskstep == 0) + { + + #ifdef STARDESTR + + #ifdef ADD_BH2 + + inp = fopen("mass-bh.con","r"); + fscanf(inp,"%lE %lE %d %lE %d", &tmp, &m_bh1, &num_bh1, &m_bh2, &num_bh2); + fclose(inp); + printf("%.8E \t %.8E %06d \t %.8E %06d \n", tmp, m_bh1, num_bh1, m_bh2, num_bh2); + fflush(stdout); + + #else + + #ifdef ADD_BH1 + + inp = fopen("mass-bh.con","r"); + fscanf(inp,"%lE %lE %d", &tmp, &m_bh1, &num_bh1); + fclose(inp); + printf("%.8E \t %.8E \t %06d \n", tmp, m_bh1, num_bh1); + fflush(stdout); + + #endif // ADD_BH1 + + #endif // ADD_BH2 + + #endif // STARDESTR + + + #ifdef STARDESTR_EXT + + inp = fopen("mass-bh.con","r"); + fscanf(inp,"%lE %lE %d", &tmp, &m_bh, &num_bh); + fclose(inp); + printf("%.8E \t %.8E \t %06d \n", tmp, m_bh, num_bh); + fflush(stdout); + + #endif // STARDESTR_EXT + + + } // if (diskstep == 0) + + + get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0); + + } /* if (myRank == rootRank) */ + + + + + #ifdef STEVOL_SSE + //SSE internal parameters (see Hurley, Pols & Tout 2000) + + value1_.neta = 0.5; //Reimers mass-loss coefficent (neta*4x10^-13; 0.5 normally) + value1_.bwind = 0.0; //Binary enhanced mass loss parameter (inactive for single) + value1_.hewind = 1.0; //Helium star mass loss factor (1.0 normally) + value1_.mxns = 3.0; //Maximum NS mass (1.8, nsflag=0; 3.0, nsflag=1) + + points_.pts1 = 0.05; //Time-step parameter in evolution phase: MS (0.05) + points_.pts2 = 0.01; //Time-step parameter in evolution phase: GB, CHeB, AGB, HeGB (0.01) + points_.pts3 = 0.02; //Time-step parameter in evolution phase: HG, HeMS (0.02) + + //value4_.sigma = 190.0; //Kick velocities - standard values... + value4_.sigma = 265.0; //Kick velocities - Hobbs++2005 values... + value4_.bhflag = 1; //bhflag > 0 allows velocity kick at BH formation + + value3_.idum = 19640916; // random number seed !!! + + //BSE internal parameters (see Hurley, Pols & Tout 2002) + + flags_.ceflag = 0; //ceflag > 0 activates spin-energy correction in common-envelope (0) #defunct# + flags_.tflag = 1; //tflag > 0 activates tidal circularisation (1) + flags_.ifflag = 0; //ifflag > 0 uses WD IFMR of HPE, 1995, MNRAS, 272, 800 (0) + flags_.nsflag = 1; //nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ, 572, 407 (1) + flags_.wdflag = 1; //wdflag > 0 uses modified-Mestel cooling for WDs (0) + + value5_.beta = 0.125; //beta is wind velocity factor: proportional to vwind**2 (1/8) + value5_.xi = 1.0; //xi is the wind accretion efficiency factor (1.0) + value5_.acc2 = 1.5; //acc2 is the Bondi-Hoyle wind accretion factor (3/2) + value5_.epsnov = 0.001; //epsnov is the fraction of accreted matter retained in nova eruption (0.001) + value5_.eddfac = 1.0; //eddfac is Eddington limit factor for mass transfer (1.0) + value5_.gamma = -1.0; //gamma is the angular momentum factor for mass lost during Roche (-1.0) + + value2_.alpha1 = 1.0; //alpha1 is the common-envelope efficiency parameter (1.0) + value2_.lambda = 0.5; //lambda is the binding energy factor for common envelope evolution (0.5) + + #endif // STEVOL_SSE + + + /* Wait to all processors to finish his works... */ + MPI_Barrier(MPI_COMM_WORLD); + + /* Broadcast all useful values to all processors... */ + MPI_Bcast(&N, 1, MPI_INT, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&eps, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&eta, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&t_end, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&dt_disk, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&dt_contr, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&dt_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&time_cur, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + + #ifdef NORM + MPI_Bcast(&m_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&r_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&v_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&t_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + #endif // NORM + + #ifdef EXTPOT + + #ifdef EXTPOT_GAL + MPI_Bcast(&m_bulge, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&a_bulge, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&b_bulge, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + + MPI_Bcast(&m_disk, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&a_disk, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&b_disk, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + + MPI_Bcast(&m_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&a_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&b_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + + #endif + + #ifdef EXTPOT_BH + MPI_Bcast(&m_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&b_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + #endif + + #ifdef EXTPOT_GAL_DEH + MPI_Bcast(&m_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&r_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&g_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + #endif + + #ifdef EXTPOT_GAL_LOG + MPI_Bcast(&v_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&r_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + + MPI_Bcast(&v2_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&r2_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + #endif + + #endif // EXTPOT + + + #ifdef STARDESTR + MPI_Bcast(&m_bh1, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&num_bh1, 1, MPI_INT, rootRank, MPI_COMM_WORLD); + + MPI_Bcast(&m_bh2, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + MPI_Bcast(&num_bh2, 1, MPI_INT, rootRank, MPI_COMM_WORLD); + #endif // STARDESTR + + + #ifdef STARDESTR_EXT + MPI_Bcast(&num_bh, 1, MPI_INT, rootRank, MPI_COMM_WORLD); + #endif // STARDESTR_EXT + + + /* Wait to all processors to finish his works... */ + MPI_Barrier(MPI_COMM_WORLD); + + + + eta_s = eta/ETA_S_CORR; + eta_bh = eta/ETA_BH_CORR; + + eps2 = SQR(eps); + + dt_min = 1.0*pow(2.0, DTMINPOWER); + dt_max = 1.0*pow(2.0, DTMAXPOWER); + + if (dt_disk == dt_contr) + dt_max = dt_contr; + else + dt_max = MIN(dt_disk, dt_contr); + + if (dt_max > 1.0) dt_max = 1.0; + + t_disk = dt_disk*(1.0+floor(time_cur/dt_disk)); + t_contr = dt_contr*(1.0+floor(time_cur/dt_contr)); + t_bh = dt_bh*(1.0+floor(time_cur/dt_bh)); + + + if (myRank == rootRank) + { + printf("t_disk = %.6E t_contr = %.6E t_bh = %.6E \n", t_disk, t_contr, t_bh); + printf("\n"); + fflush(stdout); + } /* if (myRank == rootRank) */ + + #ifdef STEVOL + dt_stevol = dt_max; + t_stevol = dt_stevol*(1.0 + floorl(time_cur/dt_stevol)); + + if (myRank == rootRank) + { + printf("t_stevol = %.6E dt_stevol = %.6E \n", t_stevol, dt_stevol); + printf("\n"); + fflush(stdout); + } /* if (myRank == rootRank) */ + #endif + + #ifdef STEVOL_SSE + dt_stevol = dt_max; + t_stevol = dt_stevol*(1.0 + floorl(time_cur/dt_stevol)); + + if (myRank == rootRank) + { + printf("t_stevol = %.6E dt_stevol = %.6E \n", t_stevol, dt_stevol); + printf("\n"); + fflush(stdout); + } /* if (myRank == rootRank) */ + #endif + + + for (i=0; i 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[i_bh1] -= pot_bh1; + pot[i_bh2] -= pot_bh2; + + for (k=0;k<3;k++) + { + a[i_bh1][k] -= a_bh1[k]; + a[i_bh2][k] -= a_bh2[k]; + + adot[i_bh1][k] -= adot_bh1[k]; + adot[i_bh2][k] -= adot_bh2[k]; + } + + // 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[i_bh1] += pot_bh1; + pot[i_bh2] += pot_bh2; + + for (k=0;k<3;k++) + { + a[i_bh1][k] += a_bh1[k]; + a[i_bh2][k] += a_bh2[k]; + + adot[i_bh1][k] += adot_bh1[k]; + adot[i_bh2][k] += adot_bh2[k]; + } + + #endif // ADD_N_BH + + + #ifdef ADD_PN_BH + + // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk + + dt_bh_tmp = dt[0]; + + tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1, + m_bh2, x_bh2, v_bh2, s_bh2, + C_NB, dt_bh_tmp, usedOrNot, + a_pn1, adot_pn1, + a_pn2, adot_pn2); + + for (k=0;k<3;k++) + { + a[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k]; + a[i_bh2][k] += a_pn2[1][k] + a_pn2[2][k] + a_pn2[3][k] + a_pn2[4][k] + a_pn2[5][k] + a_pn2[6][k]; + + adot[i_bh1][k] += adot_pn1[1][k] + adot_pn1[2][k] + adot_pn1[3][k] + adot_pn1[4][k] + adot_pn1[5][k] + adot_pn1[6][k]; + adot[i_bh2][k] += adot_pn2[1][k] + adot_pn2[2][k] + adot_pn2[3][k] + adot_pn2[4][k] + adot_pn2[5][k] + adot_pn2[6][k]; + } + + if (myRank == rootRank) + { + if (tmp_i == 505) + { + printf("PN RSDIST: %.8E \t %.8E \n", Timesteps, time_cur); + fflush(stdout); + exit(-1); + } + } + + #endif // ADD_PN_BH + + #endif // ADD_BH2 + + + + #if defined(EXTPOT) || defined(GMC) + calc_ext_grav_zero(); + #endif + + #ifdef GMC + calc_gmc_self_grav(); + calc_ext_gmc_grav(); + #endif + + + /* Wait to all processors to finish his works... */ + MPI_Barrier(MPI_COMM_WORLD); + + + + // sjuda stavim DRAG... + + #ifdef STARDISK + + for (i=0; i dt_max) dt_tmp = dt_max; + + + #ifdef STARDESTR + if (m[i] == 0.0) dt_tmp = dt_max; + #endif + + #ifdef STARDESTR_EXT + if (m[i] == 0.0) dt_tmp = dt_max; + #endif + + dt[i] = dt_tmp; + + + #ifdef DT_MIN_WARNING + if (myRank == 0) + { + if (dt[i] == dt_min) + { + printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[i], ind[i]); + fflush(stdout); + } + } + #endif + + } /* i */ + + + + + #ifdef ADD_BH2 + + /* define the min. dt over all the part. and set it also for the BH... */ + + min_dt = dt[0]; + + for (i=1; i 0.1) dt[i] = min_dt; + } /* i */ + + #endif // GMC2 + + + #ifdef GMC2222 + + /* define the max. dt for the GMC... */ + + for (i=1; i 0.1) dt[i] = dt_max; + } /* i */ + + #endif // GMC2 + + + + #ifdef ACT_DEF_LL + CreateLinkList(); + #ifdef DEBUG111 + if ( myRank == rootRank ) check_linklist(0); + #endif + #endif + + + + /* Wait to all processors to finish his works... */ + MPI_Barrier(MPI_COMM_WORLD); + + /* Scatter the "local" vectors from "global" */ + MPI_Scatter(dt, n_loc, MPI_DOUBLE, dt_loc, n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); + + /* Wait to all processors to finish his works... */ + MPI_Barrier(MPI_COMM_WORLD); + + /* load the new values for particles to the local GRAPE's */ + + nj=n_loc; + + /* load the nj particles to the G6 */ + + for (j=0; j 0) + for (int i=0; i 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; + + for (k=0;k<3;k++) + { + a_act_new[i_bh1][k] -= a_bh1[k]; + a_act_new[i_bh2][k] -= a_bh2[k]; + + adot_act_new[i_bh1][k] -= adot_bh1[k]; + adot_act_new[i_bh2][k] -= adot_bh2[k]; + } + + // 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; + + for (k=0;k<3;k++) + { + a_act_new[i_bh1][k] += a_bh1[k]; + a_act_new[i_bh2][k] += a_bh2[k]; + + adot_act_new[i_bh1][k] += adot_bh1[k]; + adot_act_new[i_bh2][k] += adot_bh2[k]; + } + + #endif // ADD_N_BH + + + #ifdef ADD_PN_BH + + // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk + + dt_bh_tmp = dt[0]; + + tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1, + m_bh2, x_bh2, v_bh2, s_bh2, + C_NB, dt_bh_tmp, usedOrNot, + a_pn1, adot_pn1, + a_pn2, adot_pn2); + + for (k=0;k<3;k++) + { + a_act_new[i_bh1][k] += a_pn1[1][k] + a_pn1[2][k] + a_pn1[3][k] + a_pn1[4][k] + a_pn1[5][k] + a_pn1[6][k]; + a_act_new[i_bh2][k] += a_pn2[1][k] + a_pn2[2][k] + a_pn2[3][k] + a_pn2[4][k] + a_pn2[5][k] + a_pn2[6][k]; + + adot_act_new[i_bh1][k] += adot_pn1[1][k] + adot_pn1[2][k] + adot_pn1[3][k] + adot_pn1[4][k] + adot_pn1[5][k] + adot_pn1[6][k]; + adot_act_new[i_bh2][k] += adot_pn2[1][k] + adot_pn2[2][k] + adot_pn2[3][k] + adot_pn2[4][k] + adot_pn2[5][k] + adot_pn2[6][k]; + } + + if (myRank == rootRank) + { + if (tmp_i == 505) + { + printf("PN RSDIST: TS = %.8E \t t = %.8E \n", Timesteps, time_cur); + fflush(stdout); + exit(-1); + } + } + + #endif // ADD_PN_BH + + #endif // ADD_BH2 + + #ifdef EXTPOT + calc_ext_grav(); + #endif + + #ifdef GMC + calc_gmc_self_grav(); + calc_ext_gmc_grav(); + #endif + + + // sjuda stavim DRAG... + + #ifdef STARDISK + + #ifdef TIMING + get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); + #endif + + calc_drag_force(); + + /* E_sd tut opjaty pljusujetsja... */ + + for (i=0; i dt_min) ) + { + power = log(dt_new)/log(2.0) - 1; + + dt_tmp = pow(2.0, (double)power); + } + + + #ifdef STARDISK + + // accretion disk criteria + // sjuda stavim izmenenija DT za schot DRAG pri priblizhenii k ploskosti Z=0 ... + + z_new_drag = x_act_new[i][2] + dt_tmp * v_act_new[i][2]; + r_new_drag2 = SQR(x_act_new[i][0]) + SQR(x_act_new[i][1]); + + if (r_new_drag2 > (R_CRIT*R_CRIT)) + { + hz = HZ; + } + else + { + hz = HZ*(sqrt(r_new_drag2)/R_CRIT); + } + + if (fabs(x_act_new[i][2]) > hz*2.2) + { + + if (z_new_drag * x_act_new[i][2] < 0.0) + { + dt_tmp *= 0.125; + } + else + { + if (fabs(z_new_drag) < hz*1.8) + { + dt_tmp *= 0.5; + } + } + + } + + #endif + + + if ( (dt_new > 2.0*dt_tmp) && + (fmod(min_t, 2.0*dt_tmp) == 0.0) && + (2.0*dt_tmp <= dt_max) ) + { + dt_tmp *= 2.0; + } + + dt_act[i] = dt_tmp; + + t_act[i] = min_t; + + pot_act[i] = pot_act_new[i]; + + for (k=0; k<3; k++) + { + x_act[i][k] = x_act_new[i][k]; + v_act[i][k] = v_act_new[i][k]; + a_act[i][k] = a_act_new[i][k]; + adot_act[i][k] = adot_act_new[i][k]; + } /* k */ + + + #ifdef DT_MIN_WARNING + if (myRank == 0) + { + if (dt_act[i] == dt_min) + { + printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]); + fflush(stdout); + } + } + #endif + + + + } /* i */ + + + + + /* define the min. dt over all the act. part. and set it also for the BH... */ + + #ifdef ADD_BH2 + + min_dt = dt_act[0]; + + for (i=1; i 0.1) dt_act[i] = min_dt; + } /* i */ + + + #endif // GMC2 + + + #ifdef GMC2222 + + /* define the max. dt for the GMC... */ + + for (i=1; i 0.1) dt_act[i] = dt_max; + } /* i */ + + #endif // GMC2 + + + + + + + #ifdef BBH_INF + if (myRank == rootRank) + { + + out = fopen("bbh.inf","a"); + + m_bh1 = m_act[i_bh1]; + m_bh2 = m_act[i_bh2]; + + for (k=0;k<3;k++) + { + x_bh1[k] = x_act[i_bh1][k]; + x_bh2[k] = x_act[i_bh2][k]; + + v_bh1[k] = v_act[i_bh1][k]; + v_bh2[k] = v_act[i_bh2][k]; + } + + for (k=0;k<3;k++) + { + x_bbhc[k] = (m_bh1*x_bh1[k] + m_bh2*x_bh2[k])/(m_bh1 + m_bh2); + v_bbhc[k] = (m_bh1*v_bh1[k] + m_bh2*v_bh2[k])/(m_bh1 + m_bh2); + } + + DR2 = SQR(x_bh1[0] - x_bh2[0]) + SQR(x_bh1[1] - x_bh2[1]) + SQR(x_bh1[2] - x_bh2[2]); + DV2 = SQR(v_bh1[0] - v_bh2[0]) + SQR(v_bh1[1] - v_bh2[1]) + SQR(v_bh1[2] - v_bh2[2]); + + EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2; + + SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB; + SEMI_a2 = SQR(SEMI_a); + + for (i=2; i= t_stevol) + { + + /* Define the current mass of all star particles... */ + + for (i=0; i= t_stevol) */ + + #ifdef TIMING + get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst); + DT_STEVOL += (CPU_tmp_user - CPU_tmp_user0); + #endif + + #endif + + + /* STEVOL_SSE routine on all the nodes */ + + #ifdef STEVOL_SSE + + + #ifdef TIMING + get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); + #endif + + if (time_cur >= t_stevol) + { + + /* Define the current mass of all star particles... */ + + for (i=0; i= t_stevol) */ + + #ifdef TIMING + get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst); + DT_STEVOL += (CPU_tmp_user - CPU_tmp_user0); + #endif + + + + + /* Update ALL the new m,r,v "j" part. to the GRAPE/GPU memory */ + + + for (j=0; j= t_bh) + { + + if (myRank == rootRank) + { + + #ifdef BH_OUT + /* Write BH data... */ + write_bh_data(); + #endif + + #ifdef BH_OUT_NB + /* Write BH NB data... */ + write_bh_nb_data(); + #endif + + #ifdef BH_OUT_NB_EXT + /* Write BH NB data... */ + write_bh_nb_data_ext(); + #endif + + } /* if (myRank == rootRank) */ + + t_bh += dt_bh; + } /* if (time_cur >= t_bh) */ + + + + + + if (time_cur >= t_contr) + { + + if (myRank == rootRank) + { + + #ifdef STARDESTR + + #ifdef ADD_BH2 + + out = fopen("mass-bh.dat","a"); + fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", + time_cur, m_bh1, num_bh1, m_bh2, num_bh2); + fclose(out); + + out = fopen("mass-bh.con","w"); + fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", + time_cur, m_bh1, num_bh1, m_bh2, num_bh2); + fclose(out); + + #else + + #ifdef ADD_BH1 + + out = fopen("mass-bh.dat","a"); + fprintf(out,"%.8E \t %.8E %06d \n", + time_cur, m_bh1, num_bh1); + fclose(out); + + out = fopen("mass-bh.con","w"); + fprintf(out,"%.8E \t %.8E %06d \n", + time_cur, m_bh1, num_bh1); + fclose(out); + + #endif // ADD_BH1 + + #endif // ADD_BH2 + + #endif // STARDESTR + + #ifdef STARDESTR_EXT + + out = fopen("phi-GRAPE.ext","w"); + fprintf(out,"%.8E \t %.8E \n", m_bh, b_bh); + fclose(out); + + out = fopen("mass-bh.dat","a"); + fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); + fclose(out); + + out = fopen("mass-bh.con","w"); + fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); + fclose(out); + + #endif // STARDESTR_EXT + + + + #ifdef ACT_DEF_LL + #ifdef DEBUG111 + if (myRank == rootRank) check_linklist((int)Timesteps); + #endif + #endif + + + energy_contr(); + + #ifdef CMCORR111 + for (i=0; i R_LIMITS_ALL) && (tmp_rv > 0.0) ) + { + + #ifdef EXTPOT + #endif + + for (k=0;k<3;k++) + { + v[i][k] *= 1.0E-01; + } /* k */ + + } /* if */ + + } /* i */ + + + + for (j=0; j= t_contr) */ + + + + + + + + + if (time_cur >= t_disk) + { + + + if (myRank == rootRank) + { + + diskstep++; + + write_snap_data(); + + #ifdef DEBUG_extra + write_snap_extra(); + #endif + + #ifdef SPH_GAS + write_data_SPH(); + #endif + + #ifdef GMC + write_snap_GMC(); + #endif + + } /* if (myRank == rootRank) */ + + #ifdef ETICS_DUMP + sprintf(out_fname, "coeffs.%06d.%02d.dat", diskstep, myRank); + grapite_dump(out_fname, 2); + #endif + + + + + #ifdef LIMITS_NEW + + for (i=0; i 900.0 ) + { + + m[i] = 0.0; + + pot[i] = 0.0; + + #ifdef EXTPOT + pot_ext[i] = 0.0; + #endif + + for (k=0;k<3;k++) + { + x[i][k] = 1000.0 + 1.0*my_rand2(); + v[i][k] = 1.0E-06*my_rand2(); + a[i][k] = 1.0E-06*my_rand2(); + adot[i][k] = 1.0E-06*my_rand2(); + } /* k */ + + t[i] = 2.0*t_end; + + dt[i] = 0.125; + + } /* if */ + + } /* i */ + + + + for (j=0; j= t_disk) */ + + + + #ifdef CPU_TIMELIMIT + if (myRank == rootRank) + { + get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst); + tmp_cpu = CPU_time_real-CPU_time_real0; + } /* if (myRank == rootRank) */ + #endif + + + #ifdef ACT_DEF_LL + ModifyLinkList(); + + #ifdef DEBUG111 + if ( (((int)Timesteps)%10000 == 0) && (myRank == rootRank) ) check_linklist((int)Timesteps); + #endif + #endif + + } /* while (time_cur < t_end) */ + + + + /* close the local GRAPE's */ + + #ifdef SAPPORO + g6_close_(&clusterid); + #else + g6_close(clusterid); + #endif + + + /* Wait to all processors to finish his works... */ + MPI_Barrier(MPI_COMM_WORLD); + + MPI_Reduce(&g6_calls, &g6_calls_sum, 1, MPI_DOUBLE, MPI_SUM, rootRank, MPI_COMM_WORLD); + + /* Wait to all processors to finish his works... */ + MPI_Barrier(MPI_COMM_WORLD); + + + if (myRank == rootRank) + { + + /* Write some output for the timestep annalize... */ + + printf("\n"); + printf("Timesteps = %.0f Total sum of integrated part. = %.0f g6_calls on all nodes = %.0f \n", Timesteps, n_act_sum, g6_calls); + printf("\n"); + printf("Real Speed = %.3f GFlops \n", 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09); + fflush(stdout); + + + #ifdef DEBUG111 + tmp_file = fopen("n_act_distr.dat","w"); + + fprintf(tmp_file,"\n"); + fprintf(tmp_file,"Total Timesteps = %.0f Total sum of integrated part. = %.0f g6_calls on all nodes = %.0f \n", Timesteps, n_act_sum, g6_calls); + fprintf(tmp_file,"\n"); + fprintf(tmp_file,"Real Speed = %.3f GFlops \n", 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09); + + for (i=1;i