/***************************************************************************** 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 *****************************************************************************/ #include "phigrape.h" Parameters parameters; // struct Parameters { // double eps = 1E-4; // double t_end = 1; // double dt_disk = 0.125; // double dt_contr = 0.125; // double dt_bh = 0.125; // double eta = 0.01; // std::string input_file_name = "data.con"; // bool dt_min_warning = true; // DT_MIN_WARNING // int live_smbh_count = 2; // ADD_BH1 or ADD_BH2 // double live_smbh_custom_eps = 0; // ADD_N_BH // bool live_smbh_output = true; // BH_OUT // bool live_smbh_neighbor_output = true; // BH_OUT_NB // NOTE needs `live_smbh_neighbor_number` // int live_smbh_neighbor_number = 10; // TODO make sure it's smaller than N? or just warn if it's not // bool binary_smbh_pn = true; // ADD_PN_BH // bool binary_smbh_influence_sphere_output = true; // BBH_INF // NOTE needs `binary_smbh_influence_radius_factor` // double binary_smbh_influence_radius_factor = 3.162277660168379497918067e+03; // R_INF // } parameters; //#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 3.162277660168379497918067e+03 // 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 #include #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]; double3() {} double3(const double x, const double y, const double z) { data[0] = x; data[1] = y; data[2] = z; } double& operator[](int i) {return data[i];} double3& operator=(const double3& a) { data[0] = a.data[0]; data[1] = a.data[1]; data[2] = a.data[2]; return *this; } double3& operator+=(const double3& a) { data[0] += a.data[0]; data[1] += a.data[1]; data[2] += a.data[2]; return *this; } double3& operator-=(const double3& a) { data[0] -= a.data[0]; data[1] -= a.data[1]; data[2] -= a.data[2]; return *this; } double3& operator/=(const double& c) { data[0] /= c; data[1] /= c; data[2] /= c; return *this; } double norm() { return sqrt(data[0]*data[0]+data[1]*data[1]+data[2]*data[2]); } operator double*() {return data;} }; double3 operator*(const double& c, const double3& a) { return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c); } double3 operator*(const double3& a, const double& c) { return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c); } double3 operator/(const double3& a, const double& c) { return double3(a.data[0]/c, a.data[1]/c, a.data[2]/c); } double3 operator+(const double3& a, const double3& b) { return double3(a.data[0]+b.data[0], a.data[1]+b.data[1], a.data[2]+b.data[2]); } double3 operator-(const double3& a, const double3& b) { return double3(a.data[0]-b.data[0], a.data[1]-b.data[1], a.data[2]-b.data[2]); } 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; double L[3]; // needed in pn_bh_spin.c // Needed for things related to BHs #include "debug.h" int ind_sort[N_MAX]; double var_sort[N_MAX]; double3 xcm, vcm, mom, xdc, vdc, a2, a3, a2dot1; char processor_name[MPI_MAX_PROCESSOR_NAME], inp_fname[30], out_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; 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; 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, // NOTE there are other "halo" variables in EXTPOT_GAL_LOG x2_ij, y2_ij, z2_ij, // NOTE this variable exist also in EXTPOT_GAL_LOG 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; // NOTE different from eps_BH // NOTE there are other m_bh and b_bh defined outside this ifdef block // NOTE this is a mess. Looks like eps_bh is never used, and the m_bh from outside the block is never used. #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 DR2, tmp_r2; double3 x_bbhc, v_bbhc; double DV2, EB, SEMI_a, SEMI_a2; // #endif // BBH_INF // #ifdef ADD_N_BH // double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3]; double3 x_bh1, x_bh2, v_bh1, v_bh2; // double pot_bh1, a_bh1[3], adot_bh1[3], // pot_bh2, a_bh2[3], adot_bh2[3]; double pot_bh1, pot_bh2; double3 a_bh1, adot_bh1, a_bh2, adot_bh2; //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 // ADD_N_BH // #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]; double3 a_pn1[7], adot_pn1[7], a_pn2[7], adot_pn2[7]; 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 // ADD_N_BH #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= 0) { tmp_a_bh = sqrt( SQR(a_bh[0]) + SQR(a_bh[1]) + SQR(a_bh[2]) ); tmp_adot_bh = sqrt( SQR(adot_bh[0]) + SQR(adot_bh[1]) + SQR(adot_bh[2]) ); if (parameters.binary_smbh_spin) { tmp_a_bh_pn0 = sqrt( SQR(a_pn[0][0]) + SQR(a_pn[0][1]) + SQR(a_pn[0][2]) ); tmp_a_bh_pn1 = sqrt( SQR(a_pn[1][0]) + SQR(a_pn[1][1]) + SQR(a_pn[1][2]) ); tmp_a_bh_pn2 = sqrt( SQR(a_pn[2][0]) + SQR(a_pn[2][1]) + SQR(a_pn[2][2]) ); tmp_a_bh_pn2_5 = sqrt( SQR(a_pn[3][0]) + SQR(a_pn[3][1]) + SQR(a_pn[3][2]) ); tmp_a_bh_pn3 = sqrt( SQR(a_pn[4][0]) + SQR(a_pn[4][1]) + SQR(a_pn[4][2]) ); tmp_a_bh_pn3_5 = sqrt( SQR(a_pn[5][0]) + SQR(a_pn[5][1]) + SQR(a_pn[5][2]) ); tmp_a_bh_spin = sqrt( SQR(a_pn[6][0]) + SQR(a_pn[6][1]) + SQR(a_pn[6][2]) ); tmp_adot_bh_pn0 = sqrt( SQR(adot_pn[0][0]) + SQR(adot_pn[0][1]) + SQR(adot_pn[0][2]) ); tmp_adot_bh_pn1 = sqrt( SQR(adot_pn[1][0]) + SQR(adot_pn[1][1]) + SQR(adot_pn[1][2]) ); tmp_adot_bh_pn2 = sqrt( SQR(adot_pn[2][0]) + SQR(adot_pn[2][1]) + SQR(adot_pn[2][2]) ); tmp_adot_bh_pn2_5 = sqrt( SQR(adot_pn[3][0]) + SQR(adot_pn[3][1]) + SQR(adot_pn[3][2]) ); tmp_adot_bh_pn3 = sqrt( SQR(adot_pn[4][0]) + SQR(adot_pn[4][1]) + SQR(adot_pn[4][2]) ); tmp_adot_bh_pn3_5 = sqrt( SQR(adot_pn[5][0]) + SQR(adot_pn[5][1]) + SQR(adot_pn[5][2]) ); tmp_adot_bh_spin = sqrt( SQR(adot_pn[6][0]) + SQR(adot_pn[6][1]) + SQR(adot_pn[6][2]) ); fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t\t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \n", time_cur, m[i], x[i][0], x[i][1], x[i][2], tmp_r, v[i][0], v[i][1], v[i][2], tmp_v, pot[i], a[i][0], a[i][1], a[i][2], tmp_a, adot[i][0], adot[i][1], adot[i][2], tmp_adot, dt[i], pot_bh, a_bh[0], a_bh[1], a_bh[2], tmp_a_bh, adot_bh[0], adot_bh[1], adot_bh[2], tmp_adot_bh, a_pn[0][0], a_pn[0][1], a_pn[0][2], tmp_a_bh_pn0, adot_pn[0][0], adot_pn[0][1], adot_pn[0][2], tmp_adot_bh_pn0, a_pn[1][0], a_pn[1][1], a_pn[1][2], tmp_a_bh_pn1, adot_pn[1][0], adot_pn[1][1], adot_pn[1][2], tmp_adot_bh_pn1, a_pn[2][0], a_pn[2][1], a_pn[2][2], tmp_a_bh_pn2, adot_pn[2][0], adot_pn[2][1], adot_pn[2][2], tmp_adot_bh_pn2, a_pn[3][0], a_pn[3][1], a_pn[3][2], tmp_a_bh_pn2_5, adot_pn[3][0], adot_pn[3][1], adot_pn[3][2], tmp_adot_bh_pn2_5, a_pn[4][0], a_pn[4][1], a_pn[4][2], tmp_a_bh_pn3, adot_pn[4][0], adot_pn[4][1], adot_pn[4][2], tmp_adot_bh_pn3, a_pn[5][0], a_pn[5][1], a_pn[5][2], tmp_a_bh_pn3_5, adot_pn[5][0], adot_pn[5][1], adot_pn[5][2], tmp_adot_bh_pn3_5, a_pn[6][0], a_pn[6][1], a_pn[6][2], tmp_a_bh_spin, adot_pn[6][0], adot_pn[6][1], adot_pn[6][2], tmp_adot_bh_spin); } else { fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \t\t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \n", time_cur, m[i], x[i][0], x[i][1], x[i][2], tmp_r, v[i][0], v[i][1], v[i][2], tmp_v, pot[i], a[i][0], a[i][1], a[i][2], tmp_a, adot[i][0], adot[i][1], adot[i][2], tmp_adot, dt[i], pot_bh, a_bh[0], a_bh[1], a_bh[2], tmp_a_bh, adot_bh[0], adot_bh[1], adot_bh[2], tmp_adot_bh); } } else { fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \n", time_cur, m[i], x[i][0], x[i][1], x[i][2], tmp_r, v[i][0], v[i][1], v[i][2], tmp_v, pot[i], a[i][0], a[i][1], a[i][2], tmp_a, adot[i][0], adot[i][1], adot[i][2], tmp_adot, dt[i]); } } fprintf(out,"\n"); fclose(out); } else if (parameters.live_smbh_count == 1) { out = fopen("bh.dat", "a"); tmp_r = sqrt( SQR(x[0][0]) + SQR(x[0][1]) + SQR(x[0][2]) ); tmp_v = sqrt( SQR(v[0][0]) + SQR(v[0][1]) + SQR(v[0][2]) ); tmp_a = sqrt( SQR(a[0][0]) + SQR(a[0][1]) + SQR(a[0][2]) ); tmp_adot = sqrt( SQR(adot[0][0]) + SQR(adot[0][1]) + SQR(adot[0][2]) ); fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \n", time_cur, m[0], x[0][0], x[0][1], x[0][2], tmp_r, v[0][0], v[0][1], v[0][2], tmp_v, pot[0], a[0][0], a[0][1], a[0][2], tmp_a, adot[0][0], adot[0][1], adot[0][2], tmp_adot, dt[0]); fprintf(out,"\n"); fclose(out); } } void write_bh_nb_data() { int i_bh, nb = parameters.live_smbh_neighbor_number; double tmp, tmp_r, tmp_v; out = fopen("bh_neighbors.dat", "a"); /* 1st BH */ i_bh = 0; for (int i_bh=0; i_bh < parameters.live_smbh_count; i_bh++) { for (i=0; i 0)) { out = fopen("bh.dat", "w"); fclose(out); } if ((parameters.live_smbh_neighbor_output) && (parameters.live_smbh_count > 0)) { out = fopen("bh_neighbors.dat", "w"); fclose(out); } if (parameters.binary_smbh_influence_sphere_output) { out = fopen("bbh_inf.dat", "w"); fclose(out); } } else { // if (diskstep == 0) } // if (diskstep == 0) get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0); } /* if (myRank == rootRank) */ /* 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 /* 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) */ for (i=0; i= 0) { m_bh1 = m[i_bh1]; m_bh2 = m[i_bh2]; for (k=0;k<3;k++) { x_bh1[k] = x[i_bh1][k]; v_bh1[k] = v[i_bh1][k]; x_bh2[k] = x[i_bh2][k]; v_bh2[k] = v[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[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, parameters.live_smbh_custom_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]; } } if (parameters.binary_smbh_pn) { // 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, (double(*)[3])a_pn1, (double(*)[3])adot_pn1, (double(*)[3])a_pn2, (double(*)[3])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); } } } } #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; if (parameters.dt_min_warning && (myRank == 0)) { if (dt[i] == dt_min) { printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[i], ind[i]); fflush(stdout); } } } /* i */ if (parameters.live_smbh_count > 0) { double min_dt = *std::min_element(dt, dt+N); for (int i=0; i 0) for (int i=0; i 0) { i_bh1 = 0; i_bh2 = 1; } #endif #ifdef TIMING get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst); DT_ACT_DEF2 += (CPU_tmp_user - CPU_tmp_user0); #endif #ifdef TIMING get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); #endif for (i=0; i= 0) { m_bh1 = m_act[i_bh1]; m_bh2 = m_act[i_bh2]; x_bh1 = x_act_new[i_bh1]; v_bh1 = v_act_new[i_bh1]; x_bh2 = x_act_new[i_bh2]; v_bh2 = v_act_new[i_bh2]; // 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; a_act_new[i_bh1] -= a_bh1; a_act_new[i_bh2] -= a_bh2; adot_act_new[i_bh1] -= adot_bh1; adot_act_new[i_bh2] -= adot_bh2; // 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, parameters.live_smbh_custom_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; a_act_new[i_bh1] += a_bh1; a_act_new[i_bh2] += a_bh2; adot_act_new[i_bh1] += adot_bh1; adot_act_new[i_bh2] += adot_bh2; } if (parameters.binary_smbh_pn) { // 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, (double(*)[3])a_pn1, (double(*)[3])adot_pn1, (double(*)[3])a_pn2, (double(*)[3])adot_pn2); a_act_new[i_bh1] += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6]; a_act_new[i_bh2] += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6]; adot_act_new[i_bh1] += adot_pn1[1] + adot_pn1[2] + adot_pn1[3] + adot_pn1[4] + adot_pn1[5] + adot_pn1[6]; adot_act_new[i_bh2] += adot_pn2[1] + adot_pn2[2] + adot_pn2[3] + adot_pn2[4] + adot_pn2[5] + adot_pn2[6]; if (myRank == rootRank) { if (tmp_i == 505) { printf("PN RSDIST: TS = %.8E \t t = %.8E \n", timesteps, time_cur); fflush(stdout); exit(-1); } } } } #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 0) && (ind_act[i] < parameters.live_smbh_count)) dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); else dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); if (dt_new < dt_min) dt_tmp = dt_min; if ((dt_new < dt_tmp) && (dt_new > 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]; x_act[i] = x_act_new[i]; v_act[i] = v_act_new[i]; a_act[i] = a_act_new[i]; adot_act[i] = adot_act_new[i]; if (parameters.dt_min_warning && (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); } } } /* i */ /* define the min. dt over all the act. part. and set it also for the BH... */ if (parameters.live_smbh_count > 0) { double min_dt = *std::min_element(dt_act, dt_act+n_act); if (parameters.live_smbh_count>=1) dt_act[i_bh1] = min_dt; if (parameters.live_smbh_count==2) dt_act[i_bh2] = min_dt; } if (parameters.binary_smbh_influence_sphere_output) { if (myRank == rootRank) { out = fopen("bbh_inf.dat", "a"); m_bh1 = m_act[i_bh1]; m_bh2 = m_act[i_bh2]; x_bh1 = x_act[i_bh1]; x_bh2 = x_act[i_bh2]; v_bh1 = v_act[i_bh1]; v_bh2 = v_act[i_bh2]; x_bbhc = (m_bh1*x_bh1 + m_bh2*x_bh2)/(m_bh1 + m_bh2); v_bbhc = (m_bh1*v_bh1 + m_bh2*v_bh2)/(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_bh) { if (myRank == rootRank) { /* Write BH data... */ if (parameters.live_smbh_output) write_bh_data(); /* Write BH NB data... */ if (parameters.live_smbh_neighbor_output) write_bh_nb_data(); } /* 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) */ #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 GRAPEs */ 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; }