/***************************************************************************** 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 NORM // Physical normalization //#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 STARDESTR // disruption of stars by BH tidal forces //#define STARDESTR_EXT // disruption of stars by external BH tidal forces //#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 LIMITS // for "0" mass particles !!! //#define R_LIMITS 1.0E+03 // 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 ACT_DEF_LL #if defined(ACT_DEF_LL) && defined(ACT_DEF_GRAPITE) #error "Contradicting preprocessor flags!" #endif #include #include #include #include #include #include #include #include /* double aaa; double aaapars[5]; extern void qwerty_(double *aaa, double *aaapars); */ #define G6_NPIPE 2048 #include "grape6.h" #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) struct double3 { double data[3]; double& operator[](int i) {return data[i];} operator double*() {return data;} }; 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, mcm, rcm_mod, vcm_mod, rcm_sum=0.0, vcm_sum=0.0, eps=0.0, eps2, a2_mod, adot2_mod, dt_tmp, dt2half, dt3over6, dt4over24, dt5over120, dtinv, dt2inv, dt3inv, a0mia1, ad04plad12, ad0plad1, 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; double3 xcm, vcm, mom, xdc, vdc, a2, a3, a2dot1; char processor_name[MPI_MAX_PROCESSOR_NAME], inp_fname[30], out_fname[30], dbg_fname[30]; /* global variables */ int N, N_star, N_bh, ind[N_MAX], name[N_MAX]; double m[N_MAX], pot[N_MAX], t[N_MAX], dt[N_MAX]; double3 x[N_MAX], v[N_MAX], a[N_MAX], adot[N_MAX]; /* local variables */ int n_loc, ind_loc[N_MAX_loc]; double m_loc[N_MAX_loc], pot_loc[N_MAX_loc], t_loc[N_MAX_loc], dt_loc[N_MAX_loc]; double3 x_loc[N_MAX_loc], v_loc[N_MAX_loc], a_loc[N_MAX_loc], adot_loc[N_MAX_loc]; /* data for active particles */ int n_act, ind_act[N_MAX]; double m_act[N_MAX], pot_act[N_MAX], t_act[N_MAX], dt_act[N_MAX], pot_act_new[N_MAX], pot_act_tmp[N_MAX], pot_act_tmp_loc[N_MAX]; double3 x_act[N_MAX], v_act[N_MAX], a_act[N_MAX], adot_act[N_MAX], x_act_new[N_MAX], v_act_new[N_MAX], a_act_new[N_MAX], adot_act_new[N_MAX], a_act_tmp[N_MAX], adot_act_tmp[N_MAX], a_act_tmp_loc[N_MAX], adot_act_tmp_loc[N_MAX]; 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; int npipe=G6_NPIPE, index_i[G6_NPIPE]; double h2_i[G6_NPIPE], p_i[G6_NPIPE]; double3 x_i[G6_NPIPE], v_i[G6_NPIPE], a_i[G6_NPIPE], jerk_i[G6_NPIPE]; int new_tunit=51, new_xunit=51; int aflag=1, jflag=1, pflag=1; double ti=0.0; double3 a2by18, a1by6, aby2; /* 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, 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 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 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 #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}; #include "pn_bh_spin.c" /* 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 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; } void read_data() { inp = fopen(inp_fname,"r"); fscanf(inp,"%d \n", &diskstep); fscanf(inp,"%d \n", &N); fscanf(inp,"%lE \n", &time_cur); for (i=0; i 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) */ 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 #ifdef EXTPOT calc_ext_grav_zero(); #endif /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); /* Energy control... */ if (myRank == rootRank) { energy_contr(); } /* if (myRank == rootRank) */ #ifdef ETICS_DUMP if (diskstep==0) { sprintf(out_fname, "coeffs.%06d.%02d.dat", 0, myRank); grapite_dump(out_fname, 2); } #endif /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); /* Scatter the "local" vectors from "global" */ MPI_Scatter(pot, n_loc, MPI_DOUBLE, pot_loc, n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); MPI_Scatter(a, 3*n_loc, MPI_DOUBLE, a_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); MPI_Scatter(adot, 3*n_loc, MPI_DOUBLE, adot_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); #ifdef ETICS_CEP // First calculate the DC grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc); // Now copy it to the global particle list memcpy(x[grapite_cep_index], xdc, 3*sizeof(double)); memcpy(v[grapite_cep_index], vdc, 3*sizeof(double)); // Now copy it to the local particle list for tha appropriate rank if (myRank == grapite_cep_index/n_loc) { memcpy(x_loc[grapite_cep_index - myRank*n_loc], xdc, 3*sizeof(double)); memcpy(v_loc[grapite_cep_index - myRank*n_loc], vdc, 3*sizeof(double)); } grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]); #endif /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); /* Define initial timestep for all particles on all nodes */ for (i=0; i dt_max) dt_tmp = dt_max; 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) 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 /* correct the active particles positions etc... on all the nodes */ #ifdef TIMING get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); #endif for (i=0; i dt_min)) { power = log(dt_new)/log(2.0) - 1; dt_tmp = pow(2.0, (double)power); } 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= 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 } /* if (myRank == rootRank) */ t_bh += dt_bh; } /* if (time_cur >= t_bh) */ if (time_cur >= t_contr) { if (myRank == rootRank) { energy_contr(); /* write cont data */ write_cont_data(); /* possible OUT for timing !!! */ #ifdef TIMING out = fopen("timing.dat","a"); 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 + DT_ACT_REDUCE; fprintf(out,"%.8E \t %.6E \t %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f \t %.3f \t %.8E %.8E %.8E \t %.8E %.8E %.8E \n", time_cur, DT_TOT, 100.0*DT_ACT_DEF1/DT_TOT, 100.0*DT_ACT_DEF2/DT_TOT, 100.0*DT_ACT_DEF3/DT_TOT, 100.0*DT_ACT_PRED/DT_TOT, 100.0*DT_ACT_GRAV/DT_TOT, 100.0*DT_EXT_GRAV/DT_TOT, 100.0*DT_GMC_GRAV/DT_TOT, 100.0*DT_GMC_GMC_GRAV/DT_TOT, 100.0*DT_EXT_GMC_GRAV/DT_TOT, 100.0*DT_ACT_CORR/DT_TOT, 100.0*DT_ACT_LOAD/DT_TOT, 100.0*DT_STEVOL/DT_TOT, 100.0*DT_STARDISK/DT_TOT, 100.0*DT_STARDESTR/DT_TOT, 100.0*DT_ACT_REDUCE/DT_TOT, CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0, Timesteps, n_act_sum, 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09); fclose(out); #endif } /* if (myRank == rootRank) */ /* possible coordinate & velocity limits for ALL particles !!! */ #ifdef ETICS_CEP // We are /inside/ a control step, so all particles must be synchronized; we can safely calculate their density centre. The acceleration and jerk currently in the memory are for the predicted position of the CEP, by calling grapite_calc_center we "correct" the position and velocity, but not the gravity at that point. // First calculate the DC grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc); // Now copy it to the global particle list memcpy(x[grapite_cep_index], xdc, 3*sizeof(double)); memcpy(v[grapite_cep_index], vdc, 3*sizeof(double)); // Now copy it to the local particle list for tha appropriate rank if (myRank == grapite_cep_index/n_loc) { memcpy(x_loc[grapite_cep_index - myRank*n_loc], xdc, 3*sizeof(double)); memcpy(v_loc[grapite_cep_index - myRank*n_loc], vdc, 3*sizeof(double)); } grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]); #endif t_contr += dt_contr; } /* if (time_cur >= t_contr) */ if (time_cur >= t_disk) { if (myRank == rootRank) { diskstep++; write_snap_data(); } /* if (myRank == rootRank) */ #ifdef ETICS_DUMP sprintf(out_fname, "coeffs.%06d.%02d.dat", diskstep, myRank); grapite_dump(out_fname, 2); #endif t_disk += dt_disk; } /* if (time_cur >= t_disk) */ } /* while (time_cur < t_end) */ /* close the local GRAPE's */ g6_close(clusterid); /* 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); } /* if (myRank == rootRank) */ /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); /* Finalize the MPI work */ MPI_Finalize(); return 0; }