yohai/phigrape
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity
phigrape/phi-GRAPE.c
Yohai Meiron 953a8286eb Multiple changed mostly focusing on active particle search 
* Added optional active particle search via GRAPite.

* Fixed ETICS_DTSCF in the Makefile.

* Disabled ETICS_CEP by default.

* Renamed and improved the initialization Python script
2020-02-26 18:54:40 -05:00

7276 lines
181 KiB
C
Raw Blame History

/*****************************************************************************
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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
/*
double aaa;
double aaapars[5];
extern void qwerty_(double *aaa, double *aaapars);
*/
#ifdef SAPPORO
#include <sapporo.h>
#include <g6lib.h>
#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 <mpi.h>
/* 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<N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2]);
// for(i=0; i<N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE %lE \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2], &tmp);
#ifndef STEVOL_SSE
// for(i=0; i<N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2]);
#endif
#ifdef STEVOL
for(i=0; i<N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %d %d \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2], &m0[i], &ZZZ[i], &t0[i], &event[i], &name[i]);
#endif
#ifdef STEVOL_SSE
for(i=0; i<N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %d %lE %lE %lE %lE %lE %lE %lE \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2], &m0[i], &ZZZ[i], &t0[i], &event[i], &SSEMass[i], &SSESpin[i], &SSERad[i], &SSELum[i], &SSETemp[i], &SSEMV[i], &SSEBV[i]);
for(i=0; i<N; i++) event_old[i] = event[i];
#endif
// for(i=0; i<N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE %d \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2], &name[i]);
// for(i=0; i<N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2]);
fclose(inp);
}
/****************************************************************************/
/****************************************************************************/
void write_snap_data()
{
sprintf(out_fname,"%06d.dat",diskstep);
out = fopen(out_fname,"w");
fprintf(out,"%06d \n", diskstep);
//#ifdef STEVOL
// fprintf(out,"%06d \n", N);
//#endif
// fprintf(out,"%06d %02d %02d %06d \n", N, N_bh, N_gmc, N_star);
fprintf(out,"%07d \n", N);
fprintf(out,"%.10E \n", time_cur);
// printf("%06d \t %06d \t %.6E \n", N, diskstep, time_cur);
// fflush(stdout);
/*
for(i=0; i<N_bh; i++)
{
fprintf(out,"%06d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \t %.6E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2], eps_BH);
}
for(i=N_bh; i<N; i++)
{
fprintf(out,"%06d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \t %.6E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2], eps);
}
*/
for(i=0; i<N; i++)
{
fprintf(out,"%07d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2]);
}
#ifndef STEVOL_SSE
// for(i=0; i<N; i++)
// {
// fprintf(out,"%06d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \n",
// ind[i],
// m[i],
// x[i][0], x[i][1], x[i][2],
// v[i][0], v[i][1], v[i][2]);
// }
#endif
#ifdef STEVOL
for(i=0; i<N; i++)
{
fprintf(out,"%06d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \t %.10E %.6E %.6E \t %01d %06d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
m0[i], ZZZ[i], t0[i], event[i], name[i]);
}
#endif
#ifdef STEVOL_SSE
for(i=0; i<N; i++)
{
fprintf(out,"%06d \t %.8E \t % .8E % .8E % .8E \t % .8E % .8E % .8E \t % .8E % .8E % .8E %01d \t %.6E %.6E %.6E %.6E %.6E % .6E % .6E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
m0[i], ZZZ[i], t0[i], event[i],
SSEMass[i], SSESpin[i], SSERad[i], SSELum[i], SSETemp[i], SSEMV[i], SSEBV[i]);
}
#endif
/*
fprintf(out,"%06d %.8E % .8E % .8E % .8E % .8E % .8E % .8E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2]);
*/
/*
fprintf(out,"%06d %.8E % .8E % .8E % .8E % .8E % .8E % .8E %08d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2], name[i]);
*/
fclose(out);
}
/****************************************************************************/
/****************************************************************************/
void write_cont_data()
{
out = fopen("data.con","w");
fprintf(out,"%06d \n", diskstep);
//#ifdef STEVOL
// fprintf(out,"%06d \n", N);
//#endif
// fprintf(out,"%06d %02d %02d %06d \n", N, N_bh, N_gmc, N_star);
fprintf(out,"%07d \n", N);
fprintf(out,"%.16E \n", time_cur);
/*
for(i=0; i<N_bh; i++)
{
fprintf(out,"%06d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \t %.6E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2], eps_BH);
}
for(i=N_bh; i<N; i++)
{
fprintf(out,"%06d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \t %.6E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2], eps);
}
*/
for(i=0; i<N; i++)
{
fprintf(out,"%07d \t %.10E \t % .10E % .10E % .10E \t % .10E % .10E % .10E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2]);
}
#ifndef STEVOL_SSE
// for(i=0; i<N; i++)
// {
// fprintf(out,"%06d \t %.16E \t % .16E % .16E % .16E \t % .16E % .16E % .16E \n",
// ind[i],
// m[i],
// x[i][0], x[i][1], x[i][2],
// v[i][0], v[i][1], v[i][2]);
// }
#endif
#ifdef STEVOL
for(i=0; i<N; i++)
{
fprintf(out,"%06d \t %.16E \t % .16E % .16E % .16E \t % .16E % .16E % .16E \t %.16E %.6E %.6E \t %01d %06d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
m0[i], ZZZ[i], t0[i], event[i], name[i]);
}
#endif
#ifdef STEVOL_SSE
for(i=0; i<N; i++)
{
fprintf(out,"%06d \t %.16E \t % .16E % .16E % .16E \t % .16E % .16E % .16E \t % .16E % .16E % .16E %01d \t %.6E %.6E %.6E %.6E %.6E % .6E % .6E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
m0[i], ZZZ[i], t0[i], event[i],
SSEMass[i], SSESpin[i], SSERad[i], SSELum[i], SSETemp[i], SSEMV[i], SSEBV[i]);
}
#endif
/*
fprintf(out,"%06d %.16E % .16E % .16E % .16E % .16E % .16E % .16E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2]);
*/
/*
fprintf(out,"%06d %.16E % .16E % .16E % .16E % .16E % .16E % .16E %08d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2], name[i]);
*/
fclose(out);
}
/****************************************************************************/
#ifdef SPH_GAS
/****************************************************************************/
void read_data_SPH()
{
inp = fopen("sph_gas.inp","r");
fscanf(inp,"%d \n", &tmp_i);
fscanf(inp,"%d \n", &N_GAS);
fscanf(inp,"%lE \n", &tmp);
// for(i=0; i<N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %d \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2], &m0[i], &ZZZ[i], &t0[i], &event[i]);
for(i=0; i<N_GAS; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE \n", &ind_gas[i], &m_gas[i], &x_gas[i][0], &x_gas[i][1], &x_gas[i][2], &v_gas[i][0], &v_gas[i][1], &v_gas[i][2]);
fclose(inp);
}
/****************************************************************************/
/****************************************************************************/
void write_data_SPH()
{
sprintf(out_fname,"%06d.gas",diskstep);
out = fopen(out_fname,"w");
fprintf(out,"%06d \n", diskstep);
fprintf(out,"%06d \n", N);
fprintf(out,"%.6E \n", time_cur);
for(i=0; i<N; i++)
{
x_star = x[i][0];
y_star = x[i][1];
z_star = x[i][2];
DEF_DEN_STAR(x_star, y_star, z_star, x_gas, m_gas, h_gas, &DEN_star);
fprintf(out,"%06d %.8E % .8E % .8E % .8E \t %.8E \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
DEN_star);
}
fclose(out);
}
/****************************************************************************/
#endif
#ifdef GMC
/****************************************************************************/
void read_data_GMC()
{
inp = fopen(inp_fname_gmc,"r");
fscanf(inp,"%d \n", &tmp_i);
fscanf(inp,"%d \n", &N_gmc);
fscanf(inp,"%lE \n", &tmp);
for(i=0; i<N_gmc; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE %lE %d \n", &ind_gmc[i], &m_gmc[i], &x_gmc[i][0], &x_gmc[i][1], &x_gmc[i][2], &v_gmc[i][0], &v_gmc[i][1], &v_gmc[i][2], &eps_gmc[i], &name_gmc[i]);
fclose(inp);
}
/****************************************************************************/
/****************************************************************************/
void write_snap_GMC()
{
sprintf(out_fname,"%06d.gmc",diskstep);
out = fopen(out_fname,"w");
fprintf(out,"%06d \n", diskstep);
fprintf(out,"%06d \n", N_gmc);
fprintf(out,"%.6E \n", time_cur);
for(i=0; i<N_gmc; i++)
{
fprintf(out,"%06d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E %08d \n",
ind_gmc[i],
m_gmc[i],
x_gmc[i][0], x_gmc[i][1], x_gmc[i][2],
v_gmc[i][0], v_gmc[i][1], v_gmc[i][2],
eps_gmc[i], name_gmc[i]);
}
fclose(out);
}
/****************************************************************************/
/****************************************************************************/
void write_cont_GMC()
{
out = fopen("gmc.con","w");
fprintf(out,"%06d \n", diskstep);
fprintf(out,"%06d \n", N_gmc);
fprintf(out,"%.9E \n", time_cur);
for(i=0; i<N_gmc; i++)
{
fprintf(out,"%06d %.9E % .9E % .9E % .9E % .9E % .9E % .9E %.9E %08d \n",
ind_gmc[i],
m_gmc[i],
x_gmc[i][0], x_gmc[i][1], x_gmc[i][2],
v_gmc[i][0], v_gmc[i][1], v_gmc[i][2],
eps_gmc[i], name_gmc[i]);
}
fclose(out);
}
/****************************************************************************/
#endif // GMC
#ifdef DEBUG_extra
/****************************************************************************/
void write_snap_extra()
{
/* sorted by dt */
for(i=0; i<N; i++)
{
ind_sort[i] = i;
var_sort[i] = dt[i];
// printf("%07d %07d \t %.8E %.8E \n", i, ind_sort[i], dt[i], var_sort[i]);
// fflush(stdout);
}
my_sort(0, (N-1), var_sort, ind_sort);
// for(j=0; j<N; j++)
// {
// i = ind_sort[j];
// printf("%07d %07d \t %.8E %.8E \n", i, ind_sort[i], dt[i], var_sort[i]);
// fflush(stdout);
// }
sprintf(out_fname,"%06d.extra_dt",diskstep);
out = fopen(out_fname,"w");
fprintf(out,"%06d \n", diskstep);
// fprintf(out,"%06d %02d %02d %06d \n", N, N_bh, N_gmc, N_star);
fprintf(out,"%07d \n", N);
fprintf(out,"%.6E \n", time_cur);
for(j=0; j<N; j++)
{
i = ind_sort[j];
fprintf(out,"%07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E \t %.6E %07d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
pot[i],
a[i][0], a[i][1], a[i][2],
adot[i][0], adot[i][1], adot[i][2],
var_sort[j], j);
}
fclose(out);
/* sorted by R */
/*
for(i=0; i<N; i++)
{
tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
ind_sort[i] = i;
var_sort[i] = tmp_r;
}
my_sort(0, (N-1), var_sort, ind_sort);
sprintf(out_fname,"%06d.extra_R",diskstep);
out = fopen(out_fname,"w");
fprintf(out,"%06d \n", diskstep);
fprintf(out,"%06d \n", N);
fprintf(out,"%.6E \n", time_cur);
m_tot = 0.0;
for(j=0; j<N; j++) m_tot += m[j];
//if(myRank == rootRank)
// {
// printf("111 %.6E \n", m_tot);
// fflush(stdout);
// }
k = 0;
m_cum = 0.0;
n_cum = 0.0;
for(j=0; j<N; j++)
{
i = ind_sort[j];
fprintf(out,"%06d %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E \t %.6E %06d \n",
// fprintf(out,"%06d %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E \t %.6E %06d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
pot[i],
pot_ext[i],
a[i][0], a[i][1], a[i][2],
adot[i][0], adot[i][1], adot[i][2],
var_sort[j], j);
#ifdef LAGR_RAD
m_cum += m[i];
n_cum += 1.0;
if( m_cum >= (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<lagr_rad_N; k++) fprintf(out,"%.6E ", lagr_rad[k]);
fprintf(out,"\n");
fclose(out);
#endif
*/
/* sorted by V */
/*
for(i=0; i<N; i++)
{
tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
ind_sort[i] = i;
var_sort[i] = tmp_v;
}
my_sort(0, (N-1), var_sort, ind_sort);
sprintf(out_fname,"%06d.extra_V",diskstep);
out = fopen(out_fname,"w");
fprintf(out,"%06d \n", diskstep);
fprintf(out,"%06d \n", N);
fprintf(out,"%.6E \n", time_cur);
for(j=0; j<N; j++)
{
i = ind_sort[j];
fprintf(out,"%06d %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E \t %.6E %06d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
pot[i],
a[i][0], a[i][1], a[i][2],
adot[i][0], adot[i][1], adot[i][2],
var_sort[j], j);
}
fclose(out);
*/
/* sorted by POT */
/*
for(i=0; i<N; i++)
{
tmp_pot = pot[i];
ind_sort[i] = i;
var_sort[i] = tmp_pot;
}
my_sort(0, (N-1), var_sort, ind_sort);
sprintf(out_fname,"%06d.extra_POT",diskstep);
out = fopen(out_fname,"w");
fprintf(out,"%06d \n", diskstep);
fprintf(out,"%06d \n", N);
fprintf(out,"%.6E \n", time_cur);
for(j=0; j<N; j++)
{
i = ind_sort[j];
fprintf(out,"%06d %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E \t %.6E %06d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
pot[i],
a[i][0], a[i][1], a[i][2],
adot[i][0], adot[i][1], adot[i][2],
var_sort[j], j);
}
fclose(out);
*/
#ifdef RMAX_OUT111
/* sorted by POT_EFF */
i_bh = N-1;
for(i=0; i<N; i++)
{
R[0] = x[i][0] - x[i_bh][0];
R[1] = x[i][1] - x[i_bh][1];
R[2] = x[i][2] - x[i_bh][2];
V[0] = v[i][0] - v[i_bh][0];
V[1] = v[i][1] - v[i_bh][1];
V[2] = v[i][2] - v[i_bh][2];
L[0] = (R[1]*V[2] - R[2]*V[1]);
L[1] = (R[2]*V[0] - R[0]*V[2]);
L[2] = (R[0]*V[1] - R[1]*V[0]);
dr2 = SQR(R[0]) + SQR(R[1]) + SQR(R[2]) + SQR(eps);
dr = sqrt(dr2);
dv2 = SQR(V[0]) + SQR(V[1]) + SQR(V[2]);
dv = sqrt(dv2);
dl2 = SQR(L[0]) + SQR(L[1]) + SQR(L[2]);
dl = sqrt(dl2);
tmp_pot = pot[i] + 0.5 * dl2/dr2;
ind_sort[i] = i;
var_sort[i] = tmp_pot;
}
my_sort(0, (N-1), var_sort, ind_sort);
sprintf(out_fname,"%06d.extra_POT_EFF",diskstep);
out = fopen(out_fname,"w");
fprintf(out,"%06d \n", diskstep);
fprintf(out,"%06d \n", N);
fprintf(out,"%.6E \n", time_cur);
for(j=0; j<N; j++)
{
i = ind_sort[j];
fprintf(out,"%06d %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E \t %.6E %06d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
pot[i],
a[i][0], a[i][1], a[i][2],
adot[i][0], adot[i][1], adot[i][2],
var_sort[j], j);
}
fclose(out);
#endif // RMAX_OUT
/* sorted by A */
/*
for(i=0; i<N; i++)
{
tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) );
ind_sort[i] = i;
var_sort[i] = tmp_a;
}
my_sort(0, (N-1), var_sort, ind_sort);
sprintf(out_fname,"%06d.extra_A",diskstep);
out = fopen(out_fname,"w");
fprintf(out,"%06d \n", diskstep);
fprintf(out,"%06d \n", N);
fprintf(out,"%.6E \n", time_cur);
for(j=0; j<N; j++)
{
i = ind_sort[j];
fprintf(out,"%06d %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E \t %.6E %06d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
pot[i],
a[i][0], a[i][1], a[i][2],
adot[i][0], adot[i][1], adot[i][2],
var_sort[j], j);
}
fclose(out);
*/
/* sorted by ADOT */
/*
for(i=0; i<N; i++)
{
tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) );
ind_sort[i] = i;
var_sort[i] = tmp_adot;
}
my_sort(0, (N-1), var_sort, ind_sort);
sprintf(out_fname,"%06d.extra_ADOT",diskstep);
out = fopen(out_fname,"w");
fprintf(out,"%06d \n", diskstep);
fprintf(out,"%06d \n", N);
fprintf(out,"%.6E \n", time_cur);
for(j=0; j<N; j++)
{
i = ind_sort[j];
fprintf(out,"%06d %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E \t %.6E %06d \n",
ind[i],
m[i],
x[i][0], x[i][1], x[i][2],
v[i][0], v[i][1], v[i][2],
pot[i],
a[i][0], a[i][1], a[i][2],
adot[i][0], adot[i][1], adot[i][2],
var_sort[j], j);
}
fclose(out);
*/
}
/****************************************************************************/
#endif
/****************************************************************************/
void write_bh_data()
{
#ifdef ADD_BH2
out = fopen("bh.dat","a");
/* 1st BH */
// i=N-2;
i=0;
tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) );
tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) );
#ifdef ADD_N_BH
tmp_a_bh = sqrt( SQR(a_bh1[0]) + SQR(a_bh1[1]) + SQR(a_bh1[2]) );
tmp_adot_bh = sqrt( SQR(adot_bh1[0]) + SQR(adot_bh1[1]) + SQR(adot_bh1[2]) );
#ifdef ADD_PN_BH
tmp_a_bh_pn0 = sqrt( SQR(a_pn1[0][0]) + SQR(a_pn1[0][1]) + SQR(a_pn1[0][2]) );
tmp_a_bh_pn1 = sqrt( SQR(a_pn1[1][0]) + SQR(a_pn1[1][1]) + SQR(a_pn1[1][2]) );
tmp_a_bh_pn2 = sqrt( SQR(a_pn1[2][0]) + SQR(a_pn1[2][1]) + SQR(a_pn1[2][2]) );
tmp_a_bh_pn2_5 = sqrt( SQR(a_pn1[3][0]) + SQR(a_pn1[3][1]) + SQR(a_pn1[3][2]) );
tmp_a_bh_pn3 = sqrt( SQR(a_pn1[4][0]) + SQR(a_pn1[4][1]) + SQR(a_pn1[4][2]) );
tmp_a_bh_pn3_5 = sqrt( SQR(a_pn1[5][0]) + SQR(a_pn1[5][1]) + SQR(a_pn1[5][2]) );
tmp_a_bh_spin = sqrt( SQR(a_pn1[6][0]) + SQR(a_pn1[6][1]) + SQR(a_pn1[6][2]) );
tmp_adot_bh_pn0 = sqrt( SQR(adot_pn1[0][0]) + SQR(adot_pn1[0][1]) + SQR(adot_pn1[0][2]) );
tmp_adot_bh_pn1 = sqrt( SQR(adot_pn1[1][0]) + SQR(adot_pn1[1][1]) + SQR(adot_pn1[1][2]) );
tmp_adot_bh_pn2 = sqrt( SQR(adot_pn1[2][0]) + SQR(adot_pn1[2][1]) + SQR(adot_pn1[2][2]) );
tmp_adot_bh_pn2_5 = sqrt( SQR(adot_pn1[3][0]) + SQR(adot_pn1[3][1]) + SQR(adot_pn1[3][2]) );
tmp_adot_bh_pn3 = sqrt( SQR(adot_pn1[4][0]) + SQR(adot_pn1[4][1]) + SQR(adot_pn1[4][2]) );
tmp_adot_bh_pn3_5 = sqrt( SQR(adot_pn1[5][0]) + SQR(adot_pn1[5][1]) + SQR(adot_pn1[5][2]) );
tmp_adot_bh_spin = sqrt( SQR(adot_pn1[6][0]) + SQR(adot_pn1[6][1]) + SQR(adot_pn1[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_bh1,
a_bh1[0], a_bh1[1], a_bh1[2], tmp_a_bh,
adot_bh1[0], adot_bh1[1], adot_bh1[2], tmp_adot_bh,
a_pn1[0][0], a_pn1[0][1], a_pn1[0][2], tmp_a_bh_pn0,
adot_pn1[0][0], adot_pn1[0][1], adot_pn1[0][2], tmp_adot_bh_pn0,
a_pn1[1][0], a_pn1[1][1], a_pn1[1][2], tmp_a_bh_pn1,
adot_pn1[1][0], adot_pn1[1][1], adot_pn1[1][2], tmp_adot_bh_pn1,
a_pn1[2][0], a_pn1[2][1], a_pn1[2][2], tmp_a_bh_pn2,
adot_pn1[2][0], adot_pn1[2][1], adot_pn1[2][2], tmp_adot_bh_pn2,
a_pn1[3][0], a_pn1[3][1], a_pn1[3][2], tmp_a_bh_pn2_5,
adot_pn1[3][0], adot_pn1[3][1], adot_pn1[3][2], tmp_adot_bh_pn2_5,
a_pn1[4][0], a_pn1[4][1], a_pn1[4][2], tmp_a_bh_pn3,
adot_pn1[4][0], adot_pn1[4][1], adot_pn1[4][2], tmp_adot_bh_pn3,
a_pn1[5][0], a_pn1[5][1], a_pn1[5][2], tmp_a_bh_pn3_5,
adot_pn1[5][0], adot_pn1[5][1], adot_pn1[5][2], tmp_adot_bh_pn3_5,
a_pn1[6][0], a_pn1[6][1], a_pn1[6][2], tmp_a_bh_spin,
adot_pn1[6][0], adot_pn1[6][1], adot_pn1[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_bh1,
a_bh1[0], a_bh1[1], a_bh1[2], tmp_a_bh,
adot_bh1[0], adot_bh1[1], adot_bh1[2], tmp_adot_bh);
#endif // ADD_PN_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]);
#endif // ADD_N_BH
/* 2nd BH */
// i=N-1;
i=1;
tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) );
tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) );
#ifdef ADD_N_BH
tmp_a_bh = sqrt( SQR(a_bh2[0]) + SQR(a_bh2[1]) + SQR(a_bh2[2]) );
tmp_adot_bh = sqrt( SQR(adot_bh2[0]) + SQR(adot_bh2[1]) + SQR(adot_bh2[2]) );
#ifdef ADD_PN_BH
tmp_a_bh_pn0 = sqrt( SQR(a_pn2[0][0]) + SQR(a_pn2[0][1]) + SQR(a_pn2[0][2]) );
tmp_a_bh_pn1 = sqrt( SQR(a_pn2[1][0]) + SQR(a_pn2[1][1]) + SQR(a_pn2[1][2]) );
tmp_a_bh_pn2 = sqrt( SQR(a_pn2[2][0]) + SQR(a_pn2[2][1]) + SQR(a_pn2[2][2]) );
tmp_a_bh_pn2_5 = sqrt( SQR(a_pn2[3][0]) + SQR(a_pn2[3][1]) + SQR(a_pn2[3][2]) );
tmp_a_bh_pn3 = sqrt( SQR(a_pn2[4][0]) + SQR(a_pn2[4][1]) + SQR(a_pn2[4][2]) );
tmp_a_bh_pn3_5 = sqrt( SQR(a_pn2[5][0]) + SQR(a_pn2[5][1]) + SQR(a_pn2[5][2]) );
tmp_a_bh_spin = sqrt( SQR(a_pn2[6][0]) + SQR(a_pn2[6][1]) + SQR(a_pn2[6][2]) );
tmp_adot_bh_pn0 = sqrt( SQR(adot_pn2[0][0]) + SQR(adot_pn2[0][1]) + SQR(adot_pn2[0][2]) );
tmp_adot_bh_pn1 = sqrt( SQR(adot_pn2[1][0]) + SQR(adot_pn2[1][1]) + SQR(adot_pn2[1][2]) );
tmp_adot_bh_pn2 = sqrt( SQR(adot_pn2[2][0]) + SQR(adot_pn2[2][1]) + SQR(adot_pn2[2][2]) );
tmp_adot_bh_pn2_5 = sqrt( SQR(adot_pn2[3][0]) + SQR(adot_pn2[3][1]) + SQR(adot_pn2[3][2]) );
tmp_adot_bh_pn3 = sqrt( SQR(adot_pn2[4][0]) + SQR(adot_pn2[4][1]) + SQR(adot_pn2[4][2]) );
tmp_adot_bh_pn3_5 = sqrt( SQR(adot_pn2[5][0]) + SQR(adot_pn2[5][1]) + SQR(adot_pn2[5][2]) );
tmp_adot_bh_spin = sqrt( SQR(adot_pn2[6][0]) + SQR(adot_pn2[6][1]) + SQR(adot_pn2[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_bh2,
a_bh2[0], a_bh2[1], a_bh2[2], tmp_a_bh,
adot_bh2[0], adot_bh2[1], adot_bh2[2], tmp_adot_bh,
a_pn2[0][0], a_pn2[0][1], a_pn2[0][2], tmp_a_bh_pn0,
adot_pn2[0][0], adot_pn2[0][1], adot_pn2[0][2], tmp_adot_bh_pn0,
a_pn2[1][0], a_pn2[1][1], a_pn2[1][2], tmp_a_bh_pn1,
adot_pn2[1][0], adot_pn2[1][1], adot_pn2[1][2], tmp_adot_bh_pn1,
a_pn2[2][0], a_pn2[2][1], a_pn2[2][2], tmp_a_bh_pn2,
adot_pn2[2][0], adot_pn2[2][1], adot_pn2[2][2], tmp_adot_bh_pn2,
a_pn2[3][0], a_pn2[3][1], a_pn2[3][2], tmp_a_bh_pn2_5,
adot_pn2[3][0], adot_pn2[3][1], adot_pn2[3][2], tmp_adot_bh_pn2_5,
a_pn2[4][0], a_pn2[4][1], a_pn2[4][2], tmp_a_bh_pn3,
adot_pn2[4][0], adot_pn2[4][1], adot_pn2[4][2], tmp_adot_bh_pn3,
a_pn2[5][0], a_pn2[5][1], a_pn2[5][2], tmp_a_bh_pn3_5,
adot_pn2[5][0], adot_pn2[5][1], adot_pn2[5][2], tmp_adot_bh_pn3_5,
a_pn2[6][0], a_pn2[6][1], a_pn2[6][2], tmp_a_bh_spin,
adot_pn2[6][0], adot_pn2[6][1], adot_pn2[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_bh2,
a_bh2[0], a_bh2[1], a_bh2[2], tmp_a_bh,
adot_bh2[0], adot_bh2[1], adot_bh2[2], tmp_adot_bh);
#endif // ADD_PN_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]);
#endif // ADD_N_BH
fprintf(out,"\n");
fclose(out);
#else
#ifdef ADD_BH1
out = fopen("bh.dat","a");
// i=N-1;
i=0;
tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) );
tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][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[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);
#endif // ADD_BH1
#endif // ADD_BH2
}
/****************************************************************************/
/****************************************************************************/
#ifdef BH_OUT_NB
void write_bh_nb_data()
{
#ifdef ADD_BH2
int i_bh, nb = 10;
double tmp, tmp_r, tmp_v;
out = fopen("bh_nb.dat","a");
/* 1st BH */
// i_bh = N-2;
i_bh = 0;
for(i=0; i<N; i++)
{
tmp_r = SQR(x[i][0] - x[i_bh][0]) + SQR(x[i][1] - x[i_bh][1]) + SQR(x[i][2] - x[i_bh][2]);
ind_sort[i] = i;
var_sort[i] = tmp_r;
}
tmp = my_select(0, (N-1), (nb-1), var_sort, ind_sort);
my_sort(0, (nb-1), var_sort, ind_sort);
fprintf(out,"%.16E \t %07d \t %.8E \t % .8E % .8E % .8E \t % .8E % .8E % .8E \t",
time_cur,
i_bh,
m[i_bh],
x[i_bh][0], x[i_bh][1], x[i_bh][2],
v[i_bh][0], v[i_bh][1], v[i_bh][2]);
for(j=1; j<nb; j++)
{
i = ind_sort[j];
tmp_r = sqrt( SQR(x[i][0] - x[i_bh][0]) + SQR(x[i][1] - x[i_bh][1]) + SQR(x[i][2] - x[i_bh][2]) );
tmp_v = sqrt( SQR(v[i][0] - v[i_bh][0]) + SQR(v[i][1] - v[i_bh][1]) + SQR(v[i][2] - v[i_bh][2]) );
fprintf(out,"%02d %07d %.8E % .8E % .8E % .8E %.8E % .8E % .8E % .8E %.8E \t",
j, i,
m[i],
x[i][0], x[i][1], x[i][2], tmp_r,
v[i][0], v[i][1], v[i][2], tmp_v);
}
fprintf(out,"\n");
/* 2nd BH */
// i_bh = N-1;
i_bh = 1;
for(i=0; i<N; i++)
{
tmp_r = SQR(x[i][0] - x[i_bh][0]) + SQR(x[i][1] - x[i_bh][1]) + SQR(x[i][2] - x[i_bh][2]);
ind_sort[i] = i;
var_sort[i] = tmp_r;
}
tmp = my_select(0, (N-1), (nb-1), var_sort, ind_sort);
my_sort(0, (nb-1), var_sort, ind_sort);
fprintf(out,"%.16E \t %07d \t %.8E \t % .8E % .8E % .8E \t % .8E % .8E % .8E \t",
time_cur,
i_bh,
m[i_bh],
x[i_bh][0], x[i_bh][1], x[i_bh][2],
v[i_bh][0], v[i_bh][1], v[i_bh][2]);
for(j=1; j<nb; j++)
{
i = ind_sort[j];
tmp_r = sqrt( SQR(x[i][0] - x[i_bh][0]) + SQR(x[i][1] - x[i_bh][1]) + SQR(x[i][2] - x[i_bh][2]) );
tmp_v = sqrt( SQR(v[i][0] - v[i_bh][0]) + SQR(v[i][1] - v[i_bh][1]) + SQR(v[i][2] - v[i_bh][2]) );
fprintf(out,"%02d %07d %.8E % .8E % .8E % .8E %.8E % .8E % .8E % .8E %.8E \t",
j, i,
m[i],
x[i][0], x[i][1], x[i][2], tmp_r,
v[i][0], v[i][1], v[i][2], tmp_v);
}
fprintf(out,"\n");
fprintf(out,"\n");
fclose(out);
#else
#ifdef ADD_BH1
int i_bh, nb = 50;
double tmp, tmp_r, tmp_v;
out = fopen("bh_nb.dat","a");
// i_bh = N-1;
i_bh = 0;
for(i=0; i<N; i++)
{
tmp_r = SQR(x[i][0] - x[i_bh][0]) + SQR(x[i][1] - x[i_bh][1]) + SQR(x[i][2] - x[i_bh][2]);
ind_sort[i] = i;
var_sort[i] = tmp_r;
}
tmp = my_select(0, (N-1), (nb-1), var_sort, ind_sort);
my_sort(0, (nb-1), var_sort, ind_sort);
fprintf(out,"%.16E \t %07d \t %.8E \t % .8E % .8E % .8E \t % .8E % .8E % .8E \t",
time_cur,
i_bh,
m[i_bh],
x[i_bh][0], x[i_bh][1], x[i_bh][2],
v[i_bh][0], v[i_bh][1], v[i_bh][2]);
for(j=1; j<nb; j++)
{
i = ind_sort[j];
tmp_r = sqrt( SQR(x[i][0] - x[i_bh][0]) + SQR(x[i][1] - x[i_bh][1]) + SQR(x[i][2] - x[i_bh][2]) );
tmp_v = sqrt( SQR(v[i][0] - v[i_bh][0]) + SQR(v[i][1] - v[i_bh][1]) + SQR(v[i][2] - v[i_bh][2]) );
fprintf(out,"%02d %07d %.8E % .8E % .8E % .8E %.8E % .8E % .8E % .8E %.8E \t",
j, i,
m[i],
x[i][0], x[i][1], x[i][2], tmp_r,
v[i][0], v[i][1], v[i][2], tmp_v);
}
fprintf(out,"\n");
fclose(out);
#endif // ADD_BH1
#endif // ADD_BH2
}
#endif
/****************************************************************************/
/****************************************************************************/
#ifdef BH_OUT_NB_EXT
void write_bh_nb_data_ext()
{
int nb = 100;
double tmp, tmp_r, tmp_v;
out = fopen("bh_nb.dat","a");
for(i=0; i<N; i++)
{
tmp_r = SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]);
ind_sort[i] = i;
var_sort[i] = tmp_r;
}
tmp = my_select(0, (N-1), (nb-1), var_sort, ind_sort);
my_sort(0, (nb-1), var_sort, ind_sort);
fprintf(out,"%.8E \t ", time_cur);
for(j=1; j<nb; j++)
{
i = ind_sort[j];
tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) );
tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) );
fprintf(out,"%03d %06d %.6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E %.6E %.6E \t",
j, i, 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]);
} /* j */
fprintf(out,"\n");
fclose(out);
}
#endif
/****************************************************************************/
/****************************************************************************/
void calc_self_grav_zero()
{
/* calc the grav for the active particles */
#ifdef SAPPORO
g6_set_ti_(&clusterid, &time_cur);
#else
g6_set_ti(clusterid, time_cur);
#endif
ni = n_act;
/* define the local phi, a, adot for these active particles */
for(i=0; i<ni; i+=npipe)
{
nn = npipe;
if(ni-i < npipe) nn = ni - i;
for(ii=0; ii<nn; ii++)
{
index_i[ii] = ind_act[i+ii];
h2_i[ii] = eps2;
for(k=0;k<3;k++)
{
x_i[ii][k] = x_act[i+ii][k]; v_i[ii][k] = v_act[i+ii][k];
} /* k */
p_i[ii] = -1.0;
for(k=0;k<3;k++)
{
a_i[ii][k] = 10.0;
jerk_i[ii][k] = 100.0;
} /* k */
} /* ii */
#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<nn; ii++)
{
if(ABS(jerk_i[ii][0]) < 1.0E-05) jerk_i[ii][0] = 1.0E-05;
if(ABS(jerk_i[ii][1]) < 1.0E-05) jerk_i[ii][1] = 1.0E-05;
if(ABS(jerk_i[ii][2]) < 1.0E-05) jerk_i[ii][2] = 1.0E-05;
if(ABS(jerk_i[ii][0]) > 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<nn; ii++)
{
pot_act_tmp[i+ii] = p_i[ii];
for(k=0;k<3;k++)
{
a_act_tmp[i+ii][k] = a_i[ii][k];
adot_act_tmp[i+ii][k] = jerk_i[ii][k];
} /* k */
} /* ii */
} /* i */
/* Store the value of the local partial force etc... */
for(i=0; i<n_act; i++)
{
pot_act_tmp_loc[i] = pot_act_tmp[i];
for(k=0;k<3;k++)
{
a_act_tmp_loc[i][k] = a_act_tmp[i][k];
adot_act_tmp_loc[i][k] = adot_act_tmp[i][k];
} /* k */
} /* i */
}
/****************************************************************************/
/****************************************************************************/
#ifdef EXTPOT_SC
void read_SC_mass()
{
inp = fopen("eff0.05.tsf10.00.tab","r");
// inp = fopen("DATAFILE","r");
for(i=0; i<M_SC; i++)
{
fscanf(inp,"%lE %lE %lE \n", &r_sc[i], &m_sc[i], &p_sc[i]);
// if(myRank == rootRank)
// {
// printf("i = %06d \t %.6E %.6E %.6E \n", i, r_sc[i], m_sc[i], p_sc[i]);
// fflush(stdout);
// } /* if(myRank == rootRank) */
}
if(myRank == rootRank)
{
i = 0; printf("i = %06d \t %.6E %.6E %.6E \n", i, r_sc[i], m_sc[i], p_sc[i]);
i = M_SC-1; printf("i = %06d \t %.6E %.6E %.6E \n", i, r_sc[i], m_sc[i], p_sc[i]);
fflush(stdout);
} /* if(myRank == rootRank) */
fclose(inp);
/* for accuracy test... */
/*
for(i=0; i<M_SC; i++)
{
m_sc[i] = 0.0;
p_sc[i] = 0.0;
}
*/
}
#endif
/****************************************************************************/
/****************************************************************************/
#ifdef EXTPOT_SC
void def_SC_mass(double r, double *M_R, double *pot_out_R)
{
int ii,ii1,ii2;
double r1=0.0,r2=0.0,m1=0.0,m2=0.0,p1=0.0,p2=0.0;
double m_tmp,p_tmp,qqq;
qqq = r*M_SC/r_sc[M_SC-1];
ii1 = qqq - 1;
ii2 = qqq + 1;
// if(myRank == rootRank)
// {
// printf("\t \t %06d %06d %.6E %.6E %.6E \n", ii1, ii2, r, M_R, pot_out_R);
// fflush(stdout);
// } /* if(myRank == rootRank) */
for(ii=ii1; ii<ii2; ii++)
{
if( (r_sc[ii] < r) && (r_sc[ii+1] > 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<ni; i++)
{
pot_ext[i] = 0.0;
#ifdef GMC
pot_gmc[i] = 0.0;
#endif
}
#ifdef EXTPOT_GALXXX
for(i=0; i<ni; i++)
{
/*
x_ij = x_act[i][0]-x_ext;
y_ij = x_act[i][1]-y_ext;
z_ij = x_act[i][2]-z_ext;
vx_ij = v_act[i][0]-vx_ext;
vy_ij = v_act[i][1]-vy_ext;
vz_ij = v_act[i][2]-vz_ext;
*/
x_ij = x_act[i][0];
y_ij = x_act[i][1];
z_ij = x_act[i][2];
vx_ij = v_act[i][0];
vy_ij = v_act[i][1];
vz_ij = v_act[i][2];
x2_ij = SQR(x_ij);
y2_ij = SQR(y_ij);
z2_ij = SQR(z_ij);
/* for BULGE */
z2_tmp = z2_ij + SQR(b_bulge);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_bulge);
r_tmp = sqrt(r2_tmp);
pot_ext[i] -= m_bulge / r_tmp;
tmp = m_bulge / (r2_tmp*r_tmp);
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij * (z_tmp + a_bulge)/z_tmp;
tmp = m_bulge / (z_tmp*r2_tmp*r2_tmp*r_tmp);
adot[i][0] += tmp * (- vx_ij*z_tmp*r2_tmp
+ 3.0*x_ij*vx_ij*x_ij*z_tmp
+ 3.0*x_ij*vy_ij*y_ij*z_tmp
+ 3.0*x_ij*vz_ij*z_ij*SQR(z_tmp + a_bulge) );
adot[i][1] += tmp * (- vy_ij*z_tmp*r2_tmp
+ 3.0*y_ij*vx_ij*x_ij*z_tmp
+ 3.0*y_ij*vy_ij*y_ij*z_tmp
+ 3.0*y_ij*vz_ij*z_ij*SQR(z_tmp + a_bulge) );
adot[i][2] += tmp * (- vz_ij*(z_tmp + a_bulge)*( x2_ij*( z2_tmp*z_tmp + a_bulge*SQR(b_bulge) ) +
y2_ij*( z2_tmp*z_tmp + a_bulge*SQR(b_bulge) ) -
( 2.0*a_bulge*(SQR(z_ij) - SQR(b_bulge))*z_tmp +
2.0*SQR(z_ij*z_ij) +
SQR(b_bulge)*SQR(z_ij) -
SQR(b_bulge)*( SQR(a_bulge) + SQR(b_bulge) ) ) )
+ 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a_bulge)
+ 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a_bulge) ) / z2_tmp ;
/* for DISK */
z2_tmp = z2_ij + SQR(b_disk);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_disk);
r_tmp = sqrt(r2_tmp);
pot_ext[i] -= m_disk / r_tmp;
tmp = m_disk / (r2_tmp*r_tmp);
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij * (z_tmp + a_disk)/z_tmp;
tmp = m_disk / (z_tmp*r2_tmp*r2_tmp*r_tmp);
adot[i][0] += tmp * (- vx_ij*z_tmp*r2_tmp
+ 3.0*x_ij*vx_ij*x_ij*z_tmp
+ 3.0*x_ij*vy_ij*y_ij*z_tmp
+ 3.0*x_ij*vz_ij*z_ij*SQR(z_tmp + a_disk) );
adot[i][1] += tmp * (- vy_ij*z_tmp*r2_tmp
+ 3.0*y_ij*vx_ij*x_ij*z_tmp
+ 3.0*y_ij*vy_ij*y_ij*z_tmp
+ 3.0*y_ij*vz_ij*z_ij*SQR(z_tmp + a_disk) );
adot[i][2] += tmp * (- vz_ij*(z_tmp + a_disk)*( x2_ij*( z2_tmp*z_tmp + a_disk*SQR(b_disk) ) +
y2_ij*( z2_tmp*z_tmp + a_disk*SQR(b_disk) ) -
( 2.0*a_disk*(SQR(z_ij) - SQR(b_disk))*z_tmp +
2.0*SQR(z_ij*z_ij) +
SQR(b_disk)*SQR(z_ij) -
SQR(b_disk)*( SQR(a_disk) + SQR(b_disk) ) ) )
+ 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a_disk)
+ 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a_disk) ) / z2_tmp ;
/* for HALO */
z2_tmp = z2_ij + SQR(b_halo);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_halo);
r_tmp = sqrt(r2_tmp);
pot_ext[i] -= m_halo / r_tmp;
tmp = m_halo / (r2_tmp*r_tmp);
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij * (z_tmp + a_halo)/z_tmp;
tmp = m_halo / (z_tmp*r2_tmp*r2_tmp*r_tmp);
adot[i][0] += tmp * (- vx_ij*z_tmp*r2_tmp
+ 3.0*x_ij*vx_ij*x_ij*z_tmp
+ 3.0*x_ij*vy_ij*y_ij*z_tmp
+ 3.0*x_ij*vz_ij*z_ij*SQR(z_tmp + a_halo) );
adot[i][1] += tmp * (- vy_ij*z_tmp*r2_tmp
+ 3.0*y_ij*vx_ij*x_ij*z_tmp
+ 3.0*y_ij*vy_ij*y_ij*z_tmp
+ 3.0*y_ij*vz_ij*z_ij*SQR(z_tmp + a_halo) );
adot[i][2] += tmp * (- vz_ij*(z_tmp + a_halo)*( x2_ij*( z2_tmp*z_tmp + a_halo*SQR(b_halo) ) +
y2_ij*( z2_tmp*z_tmp + a_halo*SQR(b_halo) ) -
( 2.0*a_halo*(SQR(z_ij) - SQR(b_halo))*z_tmp +
2.0*SQR(z_ij*z_ij) +
SQR(b_halo)*SQR(z_ij) -
SQR(b_halo)*( SQR(a_halo) + SQR(b_halo) ) ) )
+ 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a_halo)
+ 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a_halo) ) / z2_tmp ;
} /* i */
#endif
#ifdef EXTPOT_GAL
for(i=0; i<ni; i++)
{
/*
x_ij = x_act[i][0]-x_ext;
y_ij = x_act[i][1]-y_ext;
z_ij = x_act[i][2]-z_ext;
vx_ij = v_act[i][0]-vx_ext;
vy_ij = v_act[i][1]-vy_ext;
vz_ij = v_act[i][2]-vz_ext;
*/
x_ij = x_act[i][0];
y_ij = x_act[i][1];
z_ij = x_act[i][2];
vx_ij = v_act[i][0];
vy_ij = v_act[i][1];
vz_ij = v_act[i][2];
x2_ij = SQR(x_ij);
y2_ij = SQR(y_ij);
z2_ij = SQR(z_ij);
/* for BULGE */
r2 = SQR(b_bulge);
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
tmp = m_bulge / r;
pot_ext[i] -= tmp;
tmp /= r2;
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij;
adot[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
/* for DISK */
z2_tmp = z2_ij + SQR(b_disk);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_disk);
r_tmp = sqrt(r2_tmp);
pot_ext[i] -= m_disk / r_tmp;
tmp = m_disk / (r2_tmp*r_tmp);
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij * (z_tmp + a_disk)/z_tmp;
tmp = m_disk / (z_tmp*r2_tmp*r2_tmp*r_tmp);
adot[i][0] += tmp * (- vx_ij*z_tmp*r2_tmp
+ 3.0*x_ij*vx_ij*x_ij*z_tmp
+ 3.0*x_ij*vy_ij*y_ij*z_tmp
+ 3.0*x_ij*vz_ij*z_ij*SQR(z_tmp + a_disk) );
adot[i][1] += tmp * (- vy_ij*z_tmp*r2_tmp
+ 3.0*y_ij*vx_ij*x_ij*z_tmp
+ 3.0*y_ij*vy_ij*y_ij*z_tmp
+ 3.0*y_ij*vz_ij*z_ij*SQR(z_tmp + a_disk) );
adot[i][2] += tmp * (- vz_ij*(z_tmp + a_disk)*( x2_ij*( z2_tmp*z_tmp + a_disk*SQR(b_disk) ) +
y2_ij*( z2_tmp*z_tmp + a_disk*SQR(b_disk) ) -
( 2.0*a_disk*(SQR(z_ij) - SQR(b_disk))*z_tmp +
2.0*SQR(z_ij*z_ij) +
SQR(b_disk)*SQR(z_ij) -
SQR(b_disk)*( SQR(a_disk) + SQR(b_disk) ) ) )
+ 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a_disk)
+ 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a_disk) ) / z2_tmp ;
/* for HALO */
r2 = SQR(b_halo);
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
tmp = m_halo / r;
pot_ext[i] -= tmp;
tmp /= r2;
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij;
adot[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
} /* i */
#endif
#ifdef EXTPOT_GAL_DEH
for(i=0; i<ni; i++)
{
x_ij = x_act[i][0];
y_ij = x_act[i][1];
z_ij = x_act[i][2];
vx_ij = v_act[i][0];
vy_ij = v_act[i][1];
vz_ij = v_act[i][2];
tmp_r = sqrt( SQR(x_ij) + SQR(y_ij) + SQR(z_ij) );
tmp_r2 = SQR(tmp_r);
tmp_r3 = POW3(tmp_r);
dum = tmp_r/(tmp_r + r_ext);
dum2 = SQR(dum);
dum3 = POW3(dum);
dum_g = pow( dum , -g_ext );
pot_ext[i] -= ( (m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
tmp = (m_ext/tmp_r3) * dum3 * dum_g;
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij;
rv_ij = ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
tmp = ( m_ext/POW3(tmp_r+r_ext) ) * dum_g;
tmp_g = ( (r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
adot[i][0] += tmp * ( tmp_g * x_ij - vx_ij );
adot[i][1] += tmp * ( tmp_g * y_ij - vy_ij );
adot[i][2] += tmp * ( tmp_g * z_ij - vz_ij );
} /* i */
#endif
/* For simple Log potential... */
#ifdef EXTPOT_GAL_LOG
for(i=0; i<ni; i++)
{
x_ij = x_act[i][0];
y_ij = x_act[i][1];
z_ij = x_act[i][2];
vx_ij = v_act[i][0];
vy_ij = v_act[i][1];
vz_ij = v_act[i][2];
r2 = SQR(x_ij) + SQR(y_ij) + SQR(z_ij);
rv_ij = 2.0 * ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
r2_r2_halo = (r2 + r2_halo);
pot_ext[i] -= 0.5*v2_halo * log(1.0 + r2/r2_halo);
tmp = v2_halo/r2_r2_halo;
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij;
tmp /= (r2 + r2_halo);
adot[i][0] += tmp * (rv_ij * x_ij - r2_r2_halo * vx_ij);
adot[i][1] += tmp * (rv_ij * y_ij - r2_r2_halo * vy_ij);
adot[i][2] += tmp * (rv_ij * z_ij - r2_r2_halo * vz_ij);
} /* i */
#endif
/* For simple Plummer potential... */
#ifdef EXTPOT_BH
for(i=0; i<ni; i++)
{
x_ij = x_act[i][0];
y_ij = x_act[i][1];
z_ij = x_act[i][2];
vx_ij = v_act[i][0];
vy_ij = v_act[i][1];
vz_ij = v_act[i][2];
// r2 = (m_bh*SQR(b_bh) + m_act[i]*eps2)/(m_bh + m_act[i]); // for STARDISK project !!!
// r2 = 0.5*(SQR(b_bh) + eps2); // for STARDISK project !!!
r2 = SQR(b_bh);
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
tmp = m_bh / r;
pot_ext[i] -= tmp;
tmp /= r2;
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij;
adot[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
} /* i */
#endif
#ifdef EXTPOT_SC
for(i=0; i<ni; i++)
{
x_ij = x_act[i][0];
y_ij = x_act[i][1];
z_ij = x_act[i][2];
vx_ij = v_act[i][0];
vy_ij = v_act[i][1];
vz_ij = v_act[i][2];
// r2 = eps2;
r2 = 0.0;
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
M_R = 0.0;
pot_out_R = 0.0;
def_SC_mass(r, &M_R, &pot_out_R); // Mass of GAS inside < r
/* for accuracy test... */
// M_R = 0.0;
// pot_out_R = 0.0;
// if(myRank == rootRank)
// {
// printf("POT_GAS: \t \t %06d \t %.6E %.6E %.6E \n", i, r, M_R, pot_out_R);
// fflush(stdout);
// } /* if(myRank == rootRank) */
tmp = M_R / r;
pot_ext[i] -= tmp + pot_out_R;
tmp /= r2;
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij;
adot[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
} /* i */
#endif
#endif // EXTPOT
#ifdef GMC
/* For simple Plummer potential... */
for(j=0; j<N_gmc; j++)
{
for(i=0; i<ni; i++)
{
x_ij = x_act[i][0]-x_gmc[j][0];
y_ij = x_act[i][1]-x_gmc[j][1];
z_ij = x_act[i][2]-x_gmc[j][2];
vx_ij = v_act[i][0]-v_gmc[j][0];
vy_ij = v_act[i][1]-v_gmc[j][1];
vz_ij = v_act[i][2]-v_gmc[j][2];
r2 = 0.5*( eps2 + SQR(eps_gmc[j]) );
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
tmp = m_gmc[j] / r;
pot_gmc[i] -= tmp;
tmp /= r2;
a[i][0] -= tmp * x_ij;
a[i][1] -= tmp * y_ij;
a[i][2] -= tmp * z_ij;
adot[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
} /* i */
} /* j */
#endif // GMC
}
/****************************************************************************/
/****************************************************************************/
void calc_self_grav()
{
/* calc the new grav for the active particles */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
#ifdef SAPPORO
g6_set_ti_(&clusterid, &min_t);
#else
g6_set_ti(clusterid, min_t);
#endif
ni = n_act;
/* define the local phi, a, adot for these active particles */
for(i=0; i<ni; i+=npipe)
{
nn = npipe;
if(ni-i < npipe) nn = ni - i;
for(ii=0; ii<nn; ii++)
{
iii = ind_act[i+ii];
index_i[ii] = iii;
h2_i[ii] = eps2;
for(k=0;k<3;k++)
{
x_i[ii][k] = x_act_new[i+ii][k];
v_i[ii][k] = v_act_new[i+ii][k];
} /* k */
p_i[ii] = pot_act_tmp_loc[iii];
for(k=0;k<3;k++)
{
a_i[ii][k] = a_act_tmp_loc[iii][k];
jerk_i[ii][k] = adot_act_tmp_loc[iii][k];
} /* k */
} /* ii */
#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
g6_calls++;
/*
for(ii=0; ii<nn; ii++)
{
if(ABS(jerk_i[ii][0]) < 1.0E-05) jerk_i[ii][0] = 1.0E-05;
if(ABS(jerk_i[ii][1]) < 1.0E-05) jerk_i[ii][1] = 1.0E-05;
if(ABS(jerk_i[ii][2]) < 1.0E-05) jerk_i[ii][2] = 1.0E-05;
if(ABS(jerk_i[ii][0]) > 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<nn; ii++)
{
pot_act_tmp[i+ii] = p_i[ii];
for(k=0;k<3;k++)
{
a_act_tmp[i+ii][k] = a_i[ii][k];
adot_act_tmp[i+ii][k] = jerk_i[ii][k];
} /* k */
} /* ii */
} /* i */
/* Store the new value of the local partial force etc... */
for(i=0; i<n_act; i++)
{
iii = ind_act[i];
pot_act_tmp_loc[iii] = pot_act_tmp[i];
for(k=0;k<3;k++)
{
a_act_tmp_loc[iii][k] = a_act_tmp[i][k];
adot_act_tmp_loc[iii][k] = adot_act_tmp[i][k];
} /* k */
} /* i */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_ACT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
#endif
}
/****************************************************************************/
/****************************************************************************/
void calc_ext_grav()
{
#ifdef EXTPOT
/* Define the external potential for all active particles on all nodes */
ni = n_act;
for(i=0; i<ni; i++)
{
pot_act_ext[i] = 0.0;
#ifdef GMC
pot_gmc[i] = 0.0;
#endif
}
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
#ifdef EXTPOT_GALXXX
for(i=0; i<ni; i++)
{
/*
x_ij = x_act_new[i][0]-x_ext;
y_ij = x_act_new[i][1]-y_ext;
z_ij = x_act_new[i][2]-z_ext;
vx_ij = v_act_new[i][0]-vx_ext;
vy_ij = v_act_new[i][1]-vy_ext;
vz_ij = v_act_new[i][2]-vz_ext;
*/
x_ij = x_act_new[i][0];
y_ij = x_act_new[i][1];
z_ij = x_act_new[i][2];
vx_ij = v_act_new[i][0];
vy_ij = v_act_new[i][1];
vz_ij = v_act_new[i][2];
x2_ij = SQR(x_ij);
y2_ij = SQR(y_ij);
z2_ij = SQR(z_ij);
/* for BULGE */
z2_tmp = z2_ij + SQR(b_bulge);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_bulge);
r_tmp = sqrt(r2_tmp);
pot_act_ext[i] -= m_bulge / r_tmp;
tmp = m_bulge / (r2_tmp*r_tmp);
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij * (z_tmp + a_bulge)/z_tmp;
tmp = m_bulge / (z_tmp*r2_tmp*r2_tmp*r_tmp);
adot_act_new[i][0] += tmp * (- vx_ij*z_tmp*r2_tmp
+ 3.0*x_ij*vx_ij*x_ij*z_tmp
+ 3.0*x_ij*vy_ij*y_ij*z_tmp
+ 3.0*x_ij*vz_ij*z_ij*SQR(z_tmp + a_bulge) );
adot_act_new[i][1] += tmp * (- vy_ij*z_tmp*r2_tmp
+ 3.0*y_ij*vx_ij*x_ij*z_tmp
+ 3.0*y_ij*vy_ij*y_ij*z_tmp
+ 3.0*y_ij*vz_ij*z_ij*SQR(z_tmp + a_bulge) );
adot_act_new[i][2] += tmp * (- vz_ij*(z_tmp + a_bulge)*( x2_ij*( z2_tmp*z_tmp + a_bulge*SQR(b_bulge) ) +
y2_ij*( z2_tmp*z_tmp + a_bulge*SQR(b_bulge) ) -
( 2.0*a_bulge*(SQR(z_ij) - SQR(b_bulge))*z_tmp +
2.0*SQR(z_ij*z_ij) +
SQR(b_bulge)*SQR(z_ij) -
SQR(b_bulge)*( SQR(a_bulge) + SQR(b_bulge) ) ) )
+ 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a_bulge)
+ 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a_bulge) ) / z2_tmp ;
/* for DISK */
z2_tmp = z2_ij + SQR(b_disk);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_disk);
r_tmp = sqrt(r2_tmp);
pot_act_ext[i] -= m_disk / r_tmp;
tmp = m_disk / (r2_tmp*r_tmp);
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij * (z_tmp + a_disk)/z_tmp;
tmp = m_disk / (z_tmp*r2_tmp*r2_tmp*r_tmp);
adot_act_new[i][0] += tmp * (- vx_ij*z_tmp*r2_tmp
+ 3.0*x_ij*vx_ij*x_ij*z_tmp
+ 3.0*x_ij*vy_ij*y_ij*z_tmp
+ 3.0*x_ij*vz_ij*z_ij*SQR(z_tmp + a_disk) );
adot_act_new[i][1] += tmp * (- vy_ij*z_tmp*r2_tmp
+ 3.0*y_ij*vx_ij*x_ij*z_tmp
+ 3.0*y_ij*vy_ij*y_ij*z_tmp
+ 3.0*y_ij*vz_ij*z_ij*SQR(z_tmp + a_disk) );
adot_act_new[i][2] += tmp * (- vz_ij*(z_tmp + a_disk)*( x2_ij*( z2_tmp*z_tmp + a_disk*SQR(b_disk) ) +
y2_ij*( z2_tmp*z_tmp + a_disk*SQR(b_disk) ) -
( 2.0*a_disk*(SQR(z_ij) - SQR(b_disk))*z_tmp +
2.0*SQR(z_ij*z_ij) +
SQR(b_disk)*SQR(z_ij) -
SQR(b_disk)*( SQR(a_disk) + SQR(b_disk) ) ) )
+ 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a_disk)
+ 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a_disk) ) / z2_tmp ;
/* for HALO */
z2_tmp = z2_ij + SQR(b_halo);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_halo);
r_tmp = sqrt(r2_tmp);
pot_act_ext[i] -= m_halo / r_tmp;
tmp = m_halo / (r2_tmp*r_tmp);
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij * (z_tmp + a_halo)/z_tmp;
tmp = m_halo / (z_tmp*r2_tmp*r2_tmp*r_tmp);
adot_act_new[i][0] += tmp * (- vx_ij*z_tmp*r2_tmp
+ 3.0*x_ij*vx_ij*x_ij*z_tmp
+ 3.0*x_ij*vy_ij*y_ij*z_tmp
+ 3.0*x_ij*vz_ij*z_ij*SQR(z_tmp + a_halo) );
adot_act_new[i][1] += tmp * (- vy_ij*z_tmp*r2_tmp
+ 3.0*y_ij*vx_ij*x_ij*z_tmp
+ 3.0*y_ij*vy_ij*y_ij*z_tmp
+ 3.0*y_ij*vz_ij*z_ij*SQR(z_tmp + a_halo) );
adot_act_new[i][2] += tmp * (- vz_ij*(z_tmp + a_halo)*( x2_ij*( z2_tmp*z_tmp + a_halo*SQR(b_halo) ) +
y2_ij*( z2_tmp*z_tmp + a_halo*SQR(b_halo) ) -
( 2.0*a_halo*(SQR(z_ij) - SQR(b_halo))*z_tmp +
2.0*SQR(z_ij*z_ij) +
SQR(b_halo)*SQR(z_ij) -
SQR(b_halo)*( SQR(a_halo) + SQR(b_halo) ) ) )
+ 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a_halo)
+ 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a_halo) ) / z2_tmp ;
} /* i */
#endif
#ifdef EXTPOT_GAL
for(i=0; i<ni; i++)
{
/*
x_ij = x_act_new[i][0]-x_ext;
y_ij = x_act_new[i][1]-y_ext;
z_ij = x_act_new[i][2]-z_ext;
vx_ij = v_act_new[i][0]-vx_ext;
vy_ij = v_act_new[i][1]-vy_ext;
vz_ij = v_act_new[i][2]-vz_ext;
*/
x_ij = x_act_new[i][0];
y_ij = x_act_new[i][1];
z_ij = x_act_new[i][2];
vx_ij = v_act_new[i][0];
vy_ij = v_act_new[i][1];
vz_ij = v_act_new[i][2];
x2_ij = SQR(x_ij);
y2_ij = SQR(y_ij);
z2_ij = SQR(z_ij);
/* for BULGE */
r2 = SQR(b_bulge);
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
tmp = m_bulge / r;
pot_act_ext[i] -= tmp;
tmp /= r2;
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij;
adot_act_new[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot_act_new[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot_act_new[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
/* for DISK */
z2_tmp = z2_ij + SQR(b_disk);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_disk);
r_tmp = sqrt(r2_tmp);
pot_act_ext[i] -= m_disk / r_tmp;
tmp = m_disk / (r2_tmp*r_tmp);
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij * (z_tmp + a_disk)/z_tmp;
tmp = m_disk / (z_tmp*r2_tmp*r2_tmp*r_tmp);
adot_act_new[i][0] += tmp * (- vx_ij*z_tmp*r2_tmp
+ 3.0*x_ij*vx_ij*x_ij*z_tmp
+ 3.0*x_ij*vy_ij*y_ij*z_tmp
+ 3.0*x_ij*vz_ij*z_ij*SQR(z_tmp + a_disk) );
adot_act_new[i][1] += tmp * (- vy_ij*z_tmp*r2_tmp
+ 3.0*y_ij*vx_ij*x_ij*z_tmp
+ 3.0*y_ij*vy_ij*y_ij*z_tmp
+ 3.0*y_ij*vz_ij*z_ij*SQR(z_tmp + a_disk) );
adot_act_new[i][2] += tmp * (- vz_ij*(z_tmp + a_disk)*( x2_ij*( z2_tmp*z_tmp + a_disk*SQR(b_disk) ) +
y2_ij*( z2_tmp*z_tmp + a_disk*SQR(b_disk) ) -
( 2.0*a_disk*(SQR(z_ij) - SQR(b_disk))*z_tmp +
2.0*SQR(z_ij*z_ij) +
SQR(b_disk)*SQR(z_ij) -
SQR(b_disk)*( SQR(a_disk) + SQR(b_disk) ) ) )
+ 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a_disk)
+ 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a_disk) ) / z2_tmp ;
/* for HALO */
r2 = SQR(b_halo);
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
tmp = m_halo / r;
pot_act_ext[i] -= tmp;
tmp /= r2;
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij;
adot_act_new[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot_act_new[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot_act_new[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
} /* i */
#endif
#ifdef EXTPOT_GAL_DEH
for(i=0; i<ni; i++)
{
x_ij = x_act_new[i][0];
y_ij = x_act_new[i][1];
z_ij = x_act_new[i][2];
vx_ij = v_act_new[i][0];
vy_ij = v_act_new[i][1];
vz_ij = v_act_new[i][2];
tmp_r = sqrt( SQR(x_ij) + SQR(y_ij) + SQR(z_ij) );
tmp_r2 = SQR(tmp_r);
tmp_r3 = POW3(tmp_r);
dum = tmp_r/(tmp_r + r_ext);
dum2 = SQR(dum);
dum3 = POW3(dum);
dum_g = pow( dum , -g_ext );
pot_act_ext[i] -= ( (m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
tmp = (m_ext/tmp_r3) * dum3 * dum_g;
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij;
rv_ij = ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
tmp = ( m_ext/POW3(tmp_r+r_ext) ) * dum_g;
tmp_g = ( (r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
adot_act_new[i][0] += tmp * ( tmp_g * x_ij - vx_ij );
adot_act_new[i][1] += tmp * ( tmp_g * y_ij - vy_ij );
adot_act_new[i][2] += tmp * ( tmp_g * z_ij - vz_ij );
} /* i */
#endif
/* For simple Log potential... */
#ifdef EXTPOT_GAL_LOG
for(i=0; i<ni; i++)
{
x_ij = x_act_new[i][0];
y_ij = x_act_new[i][1];
z_ij = x_act_new[i][2];
vx_ij = v_act_new[i][0];
vy_ij = v_act_new[i][1];
vz_ij = v_act_new[i][2];
r2 = SQR(x_ij) + SQR(y_ij) + SQR(z_ij);
rv_ij = 2.0 * ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
r2_r2_halo = (r2 + r2_halo);
pot_act_ext[i] -= 0.5*v2_halo * log(1.0 + r2/r2_halo);
tmp = v2_halo/r2_r2_halo;
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij;
tmp /= (r2 + r2_halo);
adot_act_new[i][0] += tmp * (rv_ij * x_ij - r2_r2_halo * vx_ij);
adot_act_new[i][1] += tmp * (rv_ij * y_ij - r2_r2_halo * vy_ij);
adot_act_new[i][2] += tmp * (rv_ij * z_ij - r2_r2_halo * vz_ij);
} /* i */
#endif
/* For simple Plummer potential... */
#ifdef EXTPOT_BH
for(i=0; i<ni; i++)
{
x_ij = x_act_new[i][0];
y_ij = x_act_new[i][1];
z_ij = x_act_new[i][2];
vx_ij = v_act_new[i][0];
vy_ij = v_act_new[i][1];
vz_ij = v_act_new[i][2];
// r2 = (m_bh*SQR(b_bh) + m_act[i]*eps2)/(m_bh + m_act[i]); // for STARDISK project !!!
// r2 = 0.5*(SQR(b_bh) + eps2); // for STARDISK project !!!
r2 = SQR(b_bh);
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
tmp = m_bh / r;
pot_act_ext[i] -= tmp;
tmp /= r2;
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij;
adot_act_new[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot_act_new[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot_act_new[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
} /* i */
#endif
#ifdef EXTPOT_SC
for(i=0; i<ni; i++)
{
x_ij = x_act_new[i][0];
y_ij = x_act_new[i][1];
z_ij = x_act_new[i][2];
vx_ij = v_act_new[i][0];
vy_ij = v_act_new[i][1];
vz_ij = v_act_new[i][2];
// r2 = eps2;
r2 = 0.0;
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
M_R = 0.0;
pot_out_R = 0.0;
def_SC_mass(r, &M_R, &pot_out_R); // Mass of GAS inside < r
/* for accuracy test... */
// M_R = 0.0;
// pot_out_R = 0.0;
tmp = M_R / r;
pot_act_ext[i] -= tmp + pot_out_R;
tmp /= r2;
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij;
adot_act_new[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot_act_new[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot_act_new[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
} /* i */
#endif
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
#endif
#endif // EXTPOT
#ifdef GMC
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
/* For simple Plummer potential... */
for(j=0; j<N_gmc; j++)
{
for(i=0; i<ni; i++)
{
x_ij = x_act_new[i][0]-x_gmc[j][0];
y_ij = x_act_new[i][1]-x_gmc[j][1];
z_ij = x_act_new[i][2]-x_gmc[j][2];
vx_ij = v_act_new[i][0]-v_gmc[j][0];
vy_ij = v_act_new[i][1]-v_gmc[j][1];
vz_ij = v_act_new[i][2]-v_gmc[j][2];
r2 = 0.5*( eps2 + SQR(eps_gmc[j]) );
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
tmp = m_gmc[j] / r;
pot_gmc[i] -= tmp;
tmp /= r2;
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij;
adot_act_new[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
adot_act_new[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
adot_act_new[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
} /* i */
} /* j */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_GMC_GRAV += (CPU_tmp_user - CPU_tmp_user0);
#endif
#endif // GMC
}
/****************************************************************************/
#ifdef GMC
/****************************************************************************/
void calc_gmc_self_grav()
{
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
/* For simple Plummer potential... */
for(i=0; i<N_gmc; i++)
{
pot_gmc_gmc[i] = 0.0;
a_gmc[i][0] = 0.0;
a_gmc[i][1] = 0.0;
a_gmc[i][2] = 0.0;
}
for(j=0; j<N_gmc; j++)
{
for(i=0; i<N_gmc; i++)
{
x_ij = x_gmc[i][0]-x_gmc[j][0];
y_ij = x_gmc[i][1]-x_gmc[j][1];
z_ij = x_gmc[i][2]-x_gmc[j][2];
r2 = 0.5*( SQR(eps_gmc[i]) + SQR(eps_gmc[j]) );
r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij; r = sqrt(r2);
tmp = m_gmc[j] / r;
pot_gmc_gmc[i] -= tmp;
tmp /= r2;
a_gmc[i][0] -= tmp * x_ij;
a_gmc[i][1] -= tmp * y_ij;
a_gmc[i][2] -= tmp * z_ij;
} /* i */
} /* j */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_GMC_GMC_GRAV += (CPU_tmp_user - CPU_tmp_user0);
#endif
}
/****************************************************************************/
/****************************************************************************/
void calc_ext_gmc_grav()
{
/* Define the external potential for GMC particles on all nodes */
for(i=0; i<N_gmc; i++)
{
pot_ext_gmc[i] = 0.0;
}
#ifdef EXTPOT
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
#ifdef EXTPOT_GAL
for(i=0; i<N_gmc; i++)
{
/*
x_ij = x_gmc[i][0]-x_ext;
y_ij = x_gmc[i][1]-y_ext;
z_ij = x_gmc[i][2]-z_ext;
vx_ij = v_gmc[i][0]-vx_ext;
vy_ij = v_gmc[i][1]-vy_ext;
vz_ij = v_gmc[i][2]-vz_ext;
*/
x_ij = x_gmc[i][0];
y_ij = x_gmc[i][1];
z_ij = x_gmc[i][2];
vx_ij = v_gmc[i][0];
vy_ij = v_gmc[i][1];
vz_ij = v_gmc[i][2];
x2_ij = SQR(x_ij);
y2_ij = SQR(y_ij);
z2_ij = SQR(z_ij);
/* for BULGE */
z2_tmp = z2_ij + SQR(b_bulge);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_bulge);
r_tmp = sqrt(r2_tmp);
pot_ext_gmc[i] -= m_bulge / r_tmp;
tmp = m_bulge / (r2_tmp*r_tmp);
a_gmc[i][0] -= tmp * x_ij;
a_gmc[i][1] -= tmp * y_ij;
a_gmc[i][2] -= tmp * z_ij * (z_tmp + a_bulge)/z_tmp;
/* for DISK */
z2_tmp = z2_ij + SQR(b_disk);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_disk);
r_tmp = sqrt(r2_tmp);
pot_ext_gmc[i] -= m_disk / r_tmp;
tmp = m_disk / (r2_tmp*r_tmp);
a_gmc[i][0] -= tmp * x_ij;
a_gmc[i][1] -= tmp * y_ij;
a_gmc[i][2] -= tmp * z_ij * (z_tmp + a_disk)/z_tmp;
/* for HALO */
z2_tmp = z2_ij + SQR(b_halo);
z_tmp = sqrt(z2_tmp);
r2_tmp = x2_ij + y2_ij + SQR(z_tmp + a_halo);
r_tmp = sqrt(r2_tmp);
pot_ext_gmc[i] -= m_halo / r_tmp;
tmp = m_halo / (r2_tmp*r_tmp);
a_gmc[i][0] -= tmp * x_ij;
a_gmc[i][1] -= tmp * y_ij;
a_gmc[i][2] -= tmp * z_ij * (z_tmp + a_halo)/z_tmp;
} /* i */
#endif
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_EXT_GMC_GRAV += (CPU_tmp_user - CPU_tmp_user0);
#endif
#endif // EXTPOT
}
/****************************************************************************/
/****************************************************************************/
void GMC_kick05()
{
/* Kick(dt/2) */
for(i=0;i<N_gmc;i++)
{
v_gmc[i][0] = v_gmc[i][0] + a_gmc[i][0] * dt_gmc/2.0;
v_gmc[i][1] = v_gmc[i][1] + a_gmc[i][1] * dt_gmc/2.0;
v_gmc[i][2] = v_gmc[i][2] + a_gmc[i][2] * dt_gmc/2.0;
}
}
void GMC_kick()
{
/* Kick(dt) */
for(i=0;i<N_gmc;i++)
{
v_gmc[i][0] = v_gmc[i][0] + a_gmc[i][0] * dt_gmc;
v_gmc[i][1] = v_gmc[i][1] + a_gmc[i][1] * dt_gmc;
v_gmc[i][2] = v_gmc[i][2] + a_gmc[i][2] * dt_gmc;
}
}
void GMC_drift()
{
/* Drift(dt) */
for(i=0;i<N_gmc;i++)
{
x_gmc[i][0] = x_gmc[i][0] + v_gmc[i][0] * dt_gmc + a_gmc[i][0] * SQR(dt_gmc)/2.0;
x_gmc[i][1] = x_gmc[i][1] + v_gmc[i][1] * dt_gmc + a_gmc[i][1] * SQR(dt_gmc)/2.0;
x_gmc[i][2] = x_gmc[i][2] + v_gmc[i][2] * dt_gmc + a_gmc[i][2] * SQR(dt_gmc)/2.0;
}
}
/****************************************************************************/
#endif // GMC
/****************************************************************************/
void energy_contr()
{
E_pot = 0.0;
for(i=0; i<N; i++) E_pot += m[i]*pot[i];
E_pot *= 0.5;
E_kin = 0.0;
for(i=0; i<N; i++) E_kin += m[i]*( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
E_kin *= 0.5;
E_pot_ext = 0.0;
#ifdef EXTPOT
for(i=0; i<N; i++) E_pot_ext += m[i]*pot_ext[i];
#endif
#ifdef GMC
E_pot_gmc = 0.0;
for(i=0; i<N; i++) E_pot_gmc += m[i]*pot_gmc[i];
#endif
#ifdef GMC
E_pot_gmc_gmc = 0.0;
E_pot_ext_gmc = 0.0;
E_kin_gmc = 0.0;
for(i=0; i<N_gmc; i++) E_pot_gmc_gmc += m_gmc[i]*pot_gmc_gmc[i];
E_pot_gmc_gmc *= 0.5;
for(i=0; i<N_gmc; i++) E_kin_gmc += m_gmc[i]*( SQR(v_gmc[i][0]) + SQR(v_gmc[i][1]) + SQR(v_gmc[i][2]) );
E_kin_gmc *= 0.5;
#ifdef EXTPOT
for(i=0; i<N_gmc; i++) E_pot_ext_gmc += m_gmc[i]*pot_ext_gmc[i];
#endif
#endif
if(Timesteps == 0.0)
{
rcm_sum = 0.0;
vcm_sum = 0.0;
}
mcm = 0.0;
rcm_mod = 0.0;
vcm_mod = 0.0;
for(k=0;k<3;k++)
{
xcm[k] = 0.0;
vcm[k] = 0.0;
} /* k */
for(i=0; i<N; i++)
{
mcm += m[i];
for(k=0;k<3;k++)
{
xcm[k] += m[i] * x[i][k];
vcm[k] += m[i] * v[i][k];
} /* k */
} /* i */
for(k=0;k<3;k++)
{
xcm[k] /= mcm;
vcm[k] /= mcm;
} /* k */
rcm_mod = sqrt(xcm[0]*xcm[0]+xcm[1]*xcm[1]+xcm[2]*xcm[2]);
vcm_mod = sqrt(vcm[0]*vcm[0]+vcm[1]*vcm[1]+vcm[2]*vcm[2]);
rcm_sum += rcm_mod;
vcm_sum += vcm_mod;
for(k=0;k<3;k++) mom[k] = 0.0;
for(i=0; i<N; i++)
{
mom[0] += m[i] * ( x[i][1] * v[i][2] - x[i][2] * v[i][1] );
mom[1] += m[i] * ( x[i][2] * v[i][0] - x[i][0] * v[i][2] );
mom[2] += m[i] * ( x[i][0] * v[i][1] - x[i][1] * v[i][0] );
} /* i */
get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst);
#ifdef GMC
E_tot = E_pot + E_kin + E_pot_ext + E_pot_gmc + E_pot_gmc_gmc + E_pot_ext_gmc + E_kin_gmc;
#else
E_tot = E_pot + E_kin + E_pot_ext;
#endif
E_tot_corr = E_tot + E_corr;
E_tot_corr_sd = E_tot + E_corr + E_sd;
// if( (Timesteps == 0.0) && (time_cur == 0.0) )
if( (Timesteps == 0.0) )
{
/* initialize the E_tot_0 + etc... */
E_tot_0 = E_tot;
E_tot_corr_0 = E_tot;
E_tot_corr_sd_0 = E_tot;
// printf("1: %.6E %.6E \t % .8E % .8E \n", time_cur, Timesteps, E_tot, E_tot_0);
skip_con = 1;
}
else
{
/* read the E_tot_0 + etc... */
if( skip_con == 0 )
{
inp = fopen("contr.con","r");
fscanf(inp,"%lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE %lE",
&tmp, &tmp, &tmp, &tmp,
&tmp, &tmp, &tmp, &tmp, &tmp, &tmp, &tmp,
&E_tot_0, &tmp,
&E_tot_corr_0, &tmp,
&E_tot_corr_sd_0, &tmp,
&tmp, &tmp, &tmp, &tmp, &tmp, &tmp, &tmp, &tmp);
fclose(inp);
skip_con = 1;
}
}
DE_tot = E_tot - E_tot_0;
// printf("2: %.6E %.6E \t % .8E % .8E \n", time_cur, Timesteps, E_tot, E_tot_0);
DE_tot_corr = E_tot_corr - E_tot_corr_0;
DE_tot_corr_sd = E_tot_corr_sd - E_tot_corr_sd_0;
#if defined(STARDESTR) || defined(STARDESTR_EXT)
#ifdef STARDISK
printf("%.6E %.4E % .6E % .6E %.6E % .6E % .4E % .6E % .4E % .6E % .4E %.4E \n",
time_cur, Timesteps,
E_pot, E_pot_ext, E_kin, E_tot, DE_tot,
E_tot_corr, DE_tot_corr,
E_tot_corr_sd, DE_tot_corr_sd,
CPU_time_user-CPU_time_user0);
#else
printf("%.6E %.4E % .6E % .6E %.6E % .6E % .4E % .6E % .4E %.4E \n",
time_cur, Timesteps,
E_pot, E_pot_ext, E_kin, E_tot, DE_tot,
E_tot_corr, DE_tot_corr,
CPU_time_user-CPU_time_user0);
#endif // #ifdef STARDISK
#else
#ifdef GMC
printf("%.6E %.6E \t % .8E %.8E % .8E % .8E % .8E %.8E % .8E \t % .8E % .4E \t %.4E \n",
time_cur, Timesteps,
E_pot, E_kin, E_pot_ext, E_pot_gmc, E_pot_gmc_gmc, E_kin_gmc, E_pot_ext_gmc, E_tot, DE_tot,
CPU_time_user-CPU_time_user0);
#else
printf("%.3E %.3E % .4E %.4E % .4E % .4E % .4E %.2E\n",
time_cur, Timesteps,
E_pot, E_kin, E_pot_ext, E_tot, DE_tot,
CPU_time_user-CPU_time_user0);
#endif
// printf("%.6E %.6E \t % .8E %.8E % .8E % .8E % .8E %.8E % .8E \t % .8E % .4E \t %.4E \n",
// time_cur, Timesteps,
// E_pot, E_kin, E_pot_ext, E_pot_gmc, E_pot_gmc_gmc, E_kin_gmc, E_pot_ext_gmc, E_tot, DE_tot,
// CPU_time_user-CPU_time_user0);
#endif // #if defined(STARDESTR) || defined(STARDESTR_EXT)
/*
printf("%.4E % .4E % .4E % .4E % .6E % .4E % .4E % .4E % .4E % .4E % .6E % .4E % .6E % .4E %.4E \n",
time_cur, DE_tot_00, DE_tot_corr_sd,
E_pot, E_pot_ext, E_kin, E_tot, DE_tot,
e_pot_corr, e_pot_BH_corr, e_kin_corr, e_summ_corr,
E_tot_corr, DE_tot_corr,
E_tot_corr_sd, DE_tot_corr_sd,
CPU_time_user-CPU_time_user0);
*/
fflush(stdout);
/*
out = fopen("contr.dat","a");
fprintf(out,"%.8E %.8E %.8E %.8E \t % .8E % .8E % .8E % .8E % .8E % .8E % .8E \t % .8E % .8E \t % .8E % .8E \t % .8E % .8E \t % .8E % .8E \t % .8E % .8E % .8E \t %.8E %.8E %.8E \n",
time_cur, Timesteps, n_act_sum, g6_calls,
E_pot, E_kin, E_pot_ext, E_pot_gmc, E_pot_gmc_gmc, E_kin_gmc, E_pot_ext_gmc,
E_tot, DE_tot,
E_tot_corr, DE_tot_corr,
E_tot_corr_sd, DE_tot_corr_sd,
rcm_mod, vcm_mod,
mom[0], mom[1], mom[2],
CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0);
fclose(out);
out = fopen("contr.con","w");
fprintf(out,"%.8E %.0f %.0f %.0f \t % .8E % .8E % .8E % .8E % .8E % .8E % .8E % .8E % .8E \t % .8E % .8E \t % .8E % .8E \t % .8E % .8E \t % .8E % .8E % .8E \t %.8E %.8E %.8E \n",
time_cur, Timesteps, n_act_sum, g6_calls,
E_pot, E_kin, E_pot_ext, E_pot_gmc, E_pot_gmc_gmc, E_kin_gmc, E_pot_ext_gmc,
E_tot, DE_tot,
E_tot_corr, DE_tot_corr,
E_tot_corr_sd, DE_tot_corr_sd,
rcm_mod, vcm_mod,
mom[0], mom[1], mom[2],
CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0);
fclose(out);
*/
out = fopen("contr.dat","a");
fprintf(out,"%.8E \t %.8E %.8E %.8E \t % .8E % .8E % .8E % .8E % .8E \t % .8E % .8E \t % .8E % .8E % .8E \t %.8E %.8E %.8E \n",
time_cur, Timesteps, n_act_sum, g6_calls,
E_pot, E_kin, E_pot_ext,
E_tot, DE_tot,
rcm_mod, vcm_mod,
mom[0], mom[1], mom[2],
CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0);
fclose(out);
out = fopen("contr.con","w");
fprintf(out,"%.8E \t %.8E %.8E %.8E \t % .8E % .8E % .8E % .8E % .8E \t % .8E % .8E \t % .8E % .8E % .8E \t %.8E %.8E %.8E \n",
time_cur, Timesteps, n_act_sum, g6_calls,
E_pot, E_kin, E_pot_ext,
E_tot, DE_tot,
rcm_mod, vcm_mod,
mom[0], mom[1], mom[2],
CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0);
fclose(out);
E_tot_0 = E_tot;
// printf("3: %.6E %.6E \t % .8E % .8E \n", time_cur, Timesteps, E_tot, E_tot_0);
E_tot_corr_0 = E_tot_corr;
E_tot_corr_sd_0 = E_tot_corr_sd;
}
/****************************************************************************/
/****************************************************************************/
void cm_corr()
{
mcm = 0.0;
for(k=0;k<3;k++)
{
xcm[k] = 0.0;
vcm[k] = 0.0;
} /* k */
for(i=0; i<N; i++)
{
mcm += m[i];
for(k=0;k<3;k++)
{
xcm[k] += m[i] * x[i][k];
vcm[k] += m[i] * v[i][k];
} /* k */
} /* i */
for(k=0;k<3;k++)
{
xcm[k] /= mcm;
vcm[k] /= mcm;
} /* k */
for(i=0; i<N; i++)
{
for(k=0;k<3;k++)
{
x[i][k] -= xcm[k];
v[i][k] -= vcm[k];
} /* k */
} /* i */
}
/****************************************************************************/
/****************************************************************************/
/****************************************************************************/
/****************************************************************************/
int main(int argc, char *argv[])
{
/* INIT the rand() !!! */
srand(19640916); /* it is just my birthday :-) */
/* srand(time(NULL)); */
/* Init MPI */
MPI_Init(&argc, &argv);
/* Define the total number of processors */
MPI_Comm_size(MPI_COMM_WORLD, &n_proc);
/* Define of the Rank of each processors */
MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
/* Define the processors names */
MPI_Get_processor_name(processor_name, &name_proc);
/* Print the Rank and the names of processors */
printf("Rank of the processor %03d on %s \n", myRank, processor_name);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/*
aaa = 0.02;
qwerty_(&aaa, aaapars);
printf("aaa = %.4E \n", aaa);
for(k=0;k<4;k++) printf("aaa[%02d] = %.4E \n", k, aaapars[k]);
printf("\n");
fflush(stdout);
*/
/* read the input parameters to the rootRank */
if(myRank == rootRank)
{
inp = fopen("phi-GRAPE.cfg","r");
#ifdef GMC
fscanf(inp,"%lE %lE %lE %lE %lE %lE %s %s", &eps, &t_end, &dt_disk, &dt_contr, &dt_bh, &eta, inp_fname, inp_fname_gmc);
#else
fscanf(inp,"%lE %lE %lE %lE %lE %lE %s", &eps, &t_end, &dt_disk, &dt_contr, &dt_bh, &eta, inp_fname);
#endif
fclose(inp);
#ifdef BBH_INF
for(i=0; i<N; i++)
{
inf_event[i] = 0;
}
#endif
// dt_contr = dt_min;
/*
eps : Plummer softening parameter (can be even 0)
t_end : end time of calculation
dt_disk : interval of snapshot files output (0xxx.dat)
dt_contr : interval for the energy control output (contr.dat)
dt_bh : interval for BH output (bh.dat & bh_nb.dat)
eta : parameter for timestep determination
inp_data : name of the input file (data.inp)
*/
/* Normalization... */
#ifdef NORM
inp = fopen("phi-GRAPE.norm","r");
fscanf(inp,"%lE %lE", &m_norm, &r_norm);
fclose(inp);
m_norm *= Msol;
r_norm *= pc;
v_norm = sqrt(G*m_norm/r_norm);
t_norm = (r_norm/v_norm);
#endif // NORM
/* Read the data for EXTernal POTential */
#ifdef EXTPOT
inp = fopen("phi-GRAPE.ext","r");
#ifdef EXTPOT_GAL
fscanf(inp,"%lE %lE %lE", &m_bulge, &a_bulge, &b_bulge);
fscanf(inp,"%lE %lE %lE", &m_disk, &a_disk, &b_disk);
fscanf(inp,"%lE %lE %lE", &m_halo, &a_halo, &b_halo);
// fscanf(inp,"%lE %lE %lE", &x_ext, &y_ext, &z_ext);
// fscanf(inp,"%lE %lE %lE", &vx_ext, &vy_ext, &vz_ext);
m_bulge *= (Msol/m_norm); a_bulge *= (kpc/r_norm); b_bulge *= (kpc/r_norm);
m_disk *= (Msol/m_norm); a_disk *= (kpc/r_norm); b_disk *= (kpc/r_norm);
m_halo *= (Msol/m_norm); a_halo *= (kpc/r_norm); b_halo *= (kpc/r_norm);
#endif
#ifdef EXTPOT_BH
fscanf(inp,"%lE %lE", &m_bh, &b_bh);
// m_bh *= (Msol/m_norm); b_bh *= (kpc/r_norm);
#endif
#ifdef EXTPOT_GAL_DEH
fscanf(inp,"%lE %lE %lE", &m_ext, &r_ext, &g_ext);
// m_ext *= (Msol/m_norm); r_ext *= (pc/r_norm);
#endif
#ifdef EXTPOT_GAL_LOG
fscanf(inp,"%lE %lE", &v_halo, &r_halo);
v_halo *= (km/v_norm); r_halo *= (pc/r_norm);
v2_halo = SQR(v_halo);
r2_halo = SQR(r_halo);
#endif
fclose(inp);
#endif // EXTPOT
/* read the global data for particles to the rootRank */
read_data();
#ifdef SPH_GAS
read_data_SPH();
// printf("%d \n",N_GAS);
#endif
#ifdef GMC
read_data_GMC();
// printf("%d \n",N_gmc);
#endif
/* possible coordinate & velocity limits for ALL particles !!! */
#ifdef LIMITS_ALL_BEG
for(i=0; i<N; i++)
{
tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
// tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
tmp_rv = x[i][0]*v[i][0] + x[i][1]*v[i][1] + x[i][2]*v[i][2];
// if( (tmp_r > 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<N; i++)
{
t[i] = time_cur;
dt[i] = dt_min;
}
n_loc = N/n_proc;
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* Broadcast the values of all particles to all processors... */
MPI_Bcast(ind, N, MPI_INT, rootRank, MPI_COMM_WORLD);
MPI_Bcast(m, N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(x, 3*N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(v, 3*N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
#ifdef SPH_GAS
MPI_Bcast(&N_GAS, 1, MPI_INT, rootRank, MPI_COMM_WORLD);
MPI_Bcast(ind_gas, N_GAS, MPI_INT, rootRank, MPI_COMM_WORLD);
MPI_Bcast(m_gas, N_GAS, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(x_gas, 3*N_GAS, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(v_gas, 3*N_GAS, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
#endif
#ifdef GMC
MPI_Bcast(&N_gmc, 1, MPI_INT, rootRank, MPI_COMM_WORLD);
MPI_Bcast(ind_gmc, N_gmc, MPI_INT, rootRank, MPI_COMM_WORLD);
MPI_Bcast(m_gmc, N_gmc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(x_gmc, 3*N_gmc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(v_gmc, 3*N_gmc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(eps_gmc, N_gmc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
#endif
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
#ifdef SPH_GAS
if(myRank == rootRank)
{
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
} /* if(myRank == rootRank) */
DEF_H(x_gas, h_gas, sosed_gas);
DEF_DEN(x_gas, m_gas, h_gas, sosed_gas, DEN_gas);
if(myRank == rootRank)
{
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
printf("DEF H for SPH GAS: %.8E [sec] \n", (CPU_tmp_user - CPU_tmp_user0));
fflush(stdout);
// write the SPH GAS densities to file "data.gas"...
if(diskstep == 0)
{
out = fopen("data.gas","w");
for(i=0; i<N_GAS; i++)
{
fprintf(out,"%06d %.8E % .8E % .8E % .8E % .8E % .8E % .8E \t %.8E \n",
ind_gas[i],
m_gas[i],
x_gas[i][0], x_gas[i][1], x_gas[i][2],
v_gas[i][0], v_gas[i][1], v_gas[i][2],
DEN_gas[i]);
}
fclose(out);
} /* if(diskstep == 0) */
} /* if(myRank == rootRank) */
#endif // SPH_GAS
#ifdef STEVOL
MPI_Bcast(m0, N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(t0, N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(ZZZ, N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(event, N, MPI_INT, rootRank, MPI_COMM_WORLD);
#endif
#ifdef STEVOL_SSE
MPI_Bcast(m0, N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(t0, N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(ZZZ, N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(event, N, MPI_INT, rootRank, MPI_COMM_WORLD);
MPI_Bcast(event_old, N, MPI_INT, rootRank, MPI_COMM_WORLD);
#endif
/* Scatter the "local" vectors from "global" */
MPI_Scatter(ind, n_loc, MPI_INT, ind_loc, n_loc, MPI_INT, rootRank, MPI_COMM_WORLD);
MPI_Scatter(m, n_loc, MPI_DOUBLE, m_loc, n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Scatter(t, n_loc, MPI_DOUBLE, t_loc, n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Scatter(dt, n_loc, MPI_DOUBLE, dt_loc, n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Scatter(x, 3*n_loc, MPI_DOUBLE, x_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Scatter(v, 3*n_loc, MPI_DOUBLE, v_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* init the local GRAPE's */
#ifdef SAPPORO
#ifdef LAOHUGPUDEVMGR
setRank(&myRank);
// printf("LAOHU myRank %02d \n", myRank);
// fflush(stdout);
#endif
g6_open_(&clusterid);
npipe = g6_npipes_();
g6_set_tunit_(&new_tunit);
g6_set_xunit_(&new_xunit);
// g6_set_overflow_flag_test_mode_(aflag, jflag, pflag);
#else
//numGPU = 1;
//clusterid = myRank % numGPU;
MPI_Comm shmcomm;
MPI_Comm_split_type(MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, &shmcomm);
MPI_Comm_size(shmcomm, &numGPU);
MPI_Comm_rank(shmcomm, &clusterid);
printf("Rank of the processor %03d : Number of GPUs %01d : Cluster ID %01d \n", myRank, numGPU, clusterid);
fflush(stdout);
g6_open(clusterid);
npipe = g6_npipes();
g6_set_tunit(new_tunit);
g6_set_xunit(new_xunit);
#ifndef YEBISU
g6_set_overflow_flag_test_mode(aflag, jflag, pflag);
#endif
#endif
#ifdef AMD
g6_initialize_jp_buffer(clusterid, 10000);
#endif
#ifdef GUINNESS
g6_initialize_jp_buffer(clusterid, n_loc);
#endif
#ifdef ETICS
grapite_read_particle_tags(N, "grapite.mask", myRank, n_loc);
grapite_set_dt_exp(dt_exp);
dt_exp = time_cur; // if we don't have a binary restart then we don't remember the coefficients, and we need to calculate them now.
grapite_set_t_exp(t_exp);
#endif
/* initial load the particles to the local GRAPE's */
nj=n_loc;
/* load the nj particles to the G6 */
for(j=0; j<nj; j++)
{
for(k=0;k<3;k++)
{
a2by18[k] = 0.0; a1by6[k] = 0.0; aby2[k] = 0.0;
} /* k */
#ifdef SAPPORO
g6_set_j_particle_(&clusterid, &j, &ind_loc[j], &t_loc[j], &dt_loc[j], &m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
#else
g6_set_j_particle(clusterid, j, ind_loc[j], t_loc[j], dt_loc[j], m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
#endif
} /* j */
#ifdef AMD
g6_flush_jp_buffer(clusterid);
#endif
#ifdef ETICS
double etics_length_scale;
if (myRank == rootRank) etics_length_scale = grapite_get_length_scale(N, m, x, v); // We don't want all ranks to do it, because they need to write a file and might confuse each other
MPI_Bcast(&etics_length_scale, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
grapite_set_length_scale(etics_length_scale);
#endif
#ifdef ETICS_CEP
// First time only: get the CEP index
grapite_cep_index = grapite_get_cep_index();
// 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));
}
double zeros[3] = {0., 0., 0.}; // We haven't calculated the force yet!
grapite_update_cep(time_cur, xdc, vdc, zeros, zeros);
#endif
/* define the all particles as a active on all the processors for the first time grav calc. */
n_act = N;
for(i=0; i<n_act; i++)
{
ind_act[i] = ind[i];
iii = ind_act[i];
m_act[i] = m[iii];
for(k=0;k<3;k++)
{
x_act[i][k] = x[iii][k]; v_act[i][k] = v[iii][k];
} /* k */
t_act[i] = t[iii];
dt_act[i] = dt[iii];
} /* i */
calc_self_grav_zero();
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* Reduce the "global" vectors from "local" on all processors) */
MPI_Allreduce(pot_act_tmp, pot, n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(a_act_tmp, a, 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(adot_act_tmp, adot, 3*n_act, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
#ifdef ADD_BH2
i_bh1 = 0;
i_bh2 = 1;
#ifdef ADD_N_BH
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,
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<N; i++)
{
for(k=0; k<3; k++)
{
a_drag[i][k] = 0.0;
adot_drag[i][k] = 0.0;
}
}
E_sd = 0.0;
calc_drag_force();
/* E_sd esty po vyhodu otsjuda... */
if( (Timesteps == 0.0) ) E_sd = 0.0;
for(i=0; i<N; i++)
{
for(k=0;k<3;k++)
{
a[i][k] += a_drag[i][k];
adot[i][k] += adot_drag[i][k];
}
}
#endif // STARDISK
/* 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<N; i++)
{
a2_mod = a[i][0]*a[i][0]+a[i][1]*a[i][1]+a[i][2]*a[i][2];
adot2_mod = adot[i][0]*adot[i][0]+adot[i][1]*adot[i][1]+adot[i][2]*adot[i][2];
if(adot2_mod == 0)
dt_tmp = eta_s;
else
dt_tmp = eta_s*sqrt(a2_mod/adot2_mod);
power = log(dt_tmp)/log(2.0) - 1;
// power = log(dt_tmp)/log(2.0);
dt_tmp = pow(2.0, (double)power);
if(dt_tmp < dt_min) dt_tmp = dt_min;
if(dt_tmp > 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<N; i++)
{
if( dt[i] < min_dt ) min_dt = dt[i];
} /* i */
dt[0] = min_dt;
dt[1] = min_dt;
#else
#ifdef ADD_BH1
/* define the min. dt over all the part. and set it also for the BH... */
min_dt = dt[0];
for(i=1; i<N; i++)
{
if( dt[i] < min_dt ) min_dt = dt[i];
} /* i */
dt[0] = min_dt;
#endif // ADD_BH1
#endif // ADD_BH2
#ifdef GMC
/* define the min. dt over all the part. and set it also for the GMC... */
min_dt = dt[0];
for(i=1; i<N; i++)
{
if( dt[i] < min_dt ) min_dt = dt[i];
} /* i */
dt_gmc = min_dt;
#endif // GMC
#ifdef GMC2111
/* define the min. dt over all the part. and set it also for the GMC... */
min_dt = dt[0];
for(i=1; i<N; i++)
{
if( dt[i] < min_dt ) min_dt = dt[i];
} /* i */
for(i=1; i<N; i++)
{
if(m[i] > 0.1) dt[i] = min_dt;
} /* i */
#endif // GMC2
#ifdef GMC2222
/* define the max. dt for the GMC... */
for(i=1; i<N; i++)
{
if(m[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<nj; j++)
{
for(k=0;k<3;k++)
{
a2by18[k] = 0.0; a1by6[k] = adot_loc[j][k]*over6; aby2[k] = a_loc[j][k]*over2;
} /* k */
#ifdef SAPPORO
g6_set_j_particle_(&clusterid, &j, &ind_loc[j], &t_loc[j], &dt_loc[j], &m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
#else
g6_set_j_particle(clusterid, j, ind_loc[j], t_loc[j], dt_loc[j], m_loc[j], a2by18, a1by6, aby2, v_loc[j], x_loc[j]);
#endif
} /* j */
#ifdef AMD
g6_flush_jp_buffer(clusterid);
#endif
#ifdef STEVOL_SSE
/* Define the current mass of all star particles... */
for(i=0; i<N; i++)
{
// printf("%04d \t %.6E \t %.6E %.6E \n", i, time_cur, m0[i], m[i]);
m[i] = m_curr(i, m0[i], ZZZ[i], t0[i], time_cur);
// printf("%04d \t %.6E \t %.6E %.6E \n", i, time_cur, m0[i], m[i]);
} /* i */
/*
mcm = 0.0;
for(i=0; i<N; i++) mcm += m[i];
printf("%.4E % .6E \n", time_cur, mcm);
fflush(stdout);
*/
if(myRank == rootRank)
{
if( (diskstep == 0) && (time_cur == 0.0) )
{
// printf("%06d \t %.6E \n", diskstep, time_cur);
// fflush(stdout);
write_snap_data();
// printf("%06d \t %.6E \n", diskstep, time_cur);
// fflush(stdout);
write_cont_data();
}
} /* if(myRank == rootRank) */
#endif
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
#ifdef DEBUG_extra
write_snap_extra();
#endif
} /* if(myRank == rootRank) */
/*
if(myRank == rootRank)
{
if(505 == 505)
{
printf("SOS: %.8E \t %.8E \n", Timesteps, time_cur);
fflush(stdout);
exit(-1);
}
}
*/
/* Get the Starting time on rootRank */
if(myRank == rootRank)
{
get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0);
get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst);
tmp_cpu = CPU_time_real-CPU_time_real0;
} /* if(myRank == rootRank) */
Timesteps = 0.0;
n_act_sum = 0.0;
for(i=1; i<N+1; i++)
{
n_act_distr[i-1] = 0.0;
}
g6_calls = 0.0;
#ifdef TIMING
DT_TOT = 0.0;
DT_ACT_DEF1 = 0.0;
DT_ACT_DEF2 = 0.0;
DT_ACT_DEF3 = 0.0;
DT_ACT_PRED = 0.0;
DT_ACT_GRAV = 0.0;
DT_EXT_GRAV = 0.0;
DT_EXT_GMC_GRAV = 0.0;
DT_GMC_GMC_GRAV = 0.0;
DT_ACT_CORR = 0.0;
DT_ACT_LOAD = 0.0;
DT_STEVOL = 0.0;
DT_STARDISK = 0.0;
DT_STARDESTR = 0.0;
DT_ACT_REDUCE = 0.0;
#endif
/* The main integration loop */
#ifdef CPU_TIMELIMIT
while( (time_cur <= t_end) || (tmp_cpu <= CPU_RUNLIMIT) )
#else
while(time_cur <= t_end)
// while( Timesteps < 10000 )
#endif
{
/* Define the minimal time and the active particles on all the nodes (exclude the ZERO masses!!!) */
#ifdef ACT_DEF_LL
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
get_act_plist();
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_ACT_DEF1 += (CPU_tmp_user - CPU_tmp_user0);
#endif
#else
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
/*
min_t = t[0]+dt[0]; // min_t = 2.0*t_end;
for(i=1; i<N; i++)
{
if(m[i]!=0.0) if( ((t[i]+dt[i]) < min_t) ) min_t = t[i]+dt[i];
}
*/
// j_loc_beg = myRank*n_loc;
// j_loc_end = (myRank+1)*n_loc;
#ifdef ACT_DEF_GRAPITE
min_t_loc = grapite_get_minimum_time();
#else
min_t_loc = t[0]+dt[0];
for(j=0; j<n_loc; j++)
{
jjj = j + myRank*n_loc;
tmp = t[jjj] + dt[jjj];
if( tmp < min_t_loc ) min_t_loc = tmp;
}
#endif
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* Reduce the "global" min_t from min_t_loc "local" on all processors) */
MPI_Allreduce(&min_t_loc, &min_t, 1, MPI_DOUBLE, MPI_MIN, 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_DEF1 += (CPU_tmp_user - CPU_tmp_user0);
#endif
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
#ifdef ACT_DEF_GRAPITE
int ind_act_loc[N_MAX], n_act_loc;
grapite_active_search(min_t, ind_act_loc, &n_act_loc);
if (myRank > 0)
for(int i=0; i<n_act_loc; i++)
ind_act_loc[i] += myRank*n_loc;
int n_act_arr[256], displs[256]; // Assuming maximum of 256 processes... seems safe.
MPI_Allgather(&n_act_loc, 1, MPI_INT, n_act_arr, 1, MPI_INT, MPI_COMM_WORLD);
n_act = n_act_arr[0];
for (int i=1; i<n_proc; i++)
n_act += n_act_arr[i];
displs[0] = 0;
for (int i=1; i<n_proc; i++)
displs[i]=displs[i-1]+n_act_arr[i-1];
MPI_Allgatherv(ind_act_loc, n_act_loc, MPI_INT, ind_act, n_act_arr, displs, MPI_INT, MPI_COMM_WORLD);
#else
n_act = 0;
//#pragma omp parallel for
for(i=0; i<N; i++)
{
// if(m[i]!=0.0)
// {
if( (t[i]+dt[i]) == min_t )
{
ind_act[n_act] = i;
n_act++;
}
// }
} /* i */
#endif // ACT_DEF_GRAPITE
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_ACT_DEF2 += (CPU_tmp_user - CPU_tmp_user0);
#endif
#endif // ACT_DEF_LL
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
for(i=0; i<n_act; i++)
{
iii = ind_act[i];
m_act[i] = m[iii];
for(k=0;k<3;k++)
{
x_act[i][k] = x[iii][k];
v_act[i][k] = v[iii][k];
} /* k */
t_act[i] = t[iii];
dt_act[i] = dt[iii];
pot_act[i] = pot[iii];
#ifdef EXTPOT
pot_act_ext[i] = pot_ext[iii];
#endif
for(k=0;k<3;k++)
{
a_act[i][k] = a[iii][k];
adot_act[i][k] = adot[iii][k];
} /* k */
} /* i */
/*
// calculate the SG for the set of active particles which is deepley INSIDE the EXT BH infl. rad.
// EXPERIMENTAL feature !!!
#ifdef STARDISK
SG_CALC = 1;
summ_tmp_r = 0.0;
for(i=0; i<n_act; i++)
{
summ_tmp_r += sqrt( SQR(x_act[i][0]) + SQR(x_act[i][1]) + SQR(x_act[i][2]) );
}
aver_tmp_r = summ_tmp_r/n_act;
if(aver_tmp_r < R_INT_CUT) SG_CALC = 0; else SG_CALC = 1;
#endif
*/
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_ACT_DEF3 += (CPU_tmp_user - CPU_tmp_user0);
#endif
#ifdef DEBUG111
if(myRank == rootRank)
{
printf("Timesteps = %.8E \t time_cur = %.8E \t min_t = %.8E \t n_act = %07d \n", Timesteps, time_cur, min_t, n_act);
fflush(stdout);
}
#endif
/* predict 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<n_act; i++)
{
dt_tmp = min_t - t_act[i];
dt2half = over2*dt_tmp*dt_tmp;
dt3over6 = over3*dt_tmp*dt2half;
for(k=0;k<3;k++)
{
x_act_new[i][k] = x_act[i][k] + v_act[i][k]*dt_tmp + a_act[i][k]*dt2half + adot_act[i][k]*dt3over6;
v_act_new[i][k] = v_act[i][k] + a_act[i][k]*dt_tmp + adot_act[i][k]*dt2half;
} /* k */
} /* i */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_ACT_PRED += (CPU_tmp_user - CPU_tmp_user0);
#endif
#ifdef GMC
GMC_drift();
GMC_kick();
#endif
/*
// calculate the SG for the set of active particles which is deepley INSIDE the EXT BH infl. rad.
// EXPERIMENTAL feature !!!
#ifdef STARDISK
if(SG_CALC == 1) calc_self_grav();
#else
calc_self_grav();
#endif
*/
calc_self_grav();
/* Reduce the "global" vectors from "local" on all the nodes */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
// if(min_t >= 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
i_bh1 = 0;
i_bh2 = 1;
#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<n_act; i++)
{
for(k=0;k<3;k++)
{
a_act_new[i][k] += a_drag[ind_act[i]][k];
adot_act_new[i][k] += adot_drag[ind_act[i]][k];
}
}
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_STARDISK += (CPU_tmp_user - CPU_tmp_user0);
#endif
#endif // STARDISK
/* 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<n_act; i++)
{
dt_tmp = min_t - t_act[i];
dt3over6 = dt_tmp*dt_tmp*dt_tmp/6.0;
dt4over24 = dt3over6*dt_tmp/4.0;
dt5over120 = dt4over24*dt_tmp/5.0;
dtinv = 1.0/dt_tmp;
dt2inv = dtinv*dtinv;
dt3inv = dt2inv*dtinv;
for(k=0; k<3; k++)
{
a0mia1 = a_act[i][k]-a_act_new[i][k];
ad04plad12 = 4.0*adot_act[i][k] + 2.0*adot_act_new[i][k];
ad0plad1 = adot_act[i][k] + adot_act_new[i][k];
a2[k] = -6.0*a0mia1*dt2inv - ad04plad12*dtinv;
a3[k] = 12.0*a0mia1*dt3inv + 6.0*ad0plad1*dt2inv;
x_act_new[i][k] += dt4over24*a2[k] + dt5over120*a3[k];
v_act_new[i][k] += dt3over6*a2[k] + dt4over24*a3[k];
} /* k */
a1abs = sqrt(a_act_new[i][0]*a_act_new[i][0]+a_act_new[i][1]*a_act_new[i][1]+a_act_new[i][2]*a_act_new[i][2]);
adot1abs = sqrt(adot_act_new[i][0]*adot_act_new[i][0]+adot_act_new[i][1]*adot_act_new[i][1]+adot_act_new[i][2]*adot_act_new[i][2]);
for(k=0; k<3; k++) a2dot1[k] = a2[k] + dt_tmp*a3[k];
a2dot1abs = sqrt(a2dot1[0]*a2dot1[0]+a2dot1[1]*a2dot1[1]+a2dot1[2]*a2dot1[2]);
a3dot1abs = sqrt(a3[0]*a3[0]+a3[1]*a3[1]+a3[2]*a3[2]);
// dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
#ifdef ADD_BH2
if( (ind_act[i] == 0) || (ind_act[i] == 1) )
{
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));
}
#else
#ifdef ADD_BH1
if( (ind_act[i] == 0) )
{
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));
}
#endif // ADD_BH1
dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
#endif // ADD_BH2
#ifdef STARDESTR
if(m_act[i] == 0.0) dt_new = dt_max;
#endif
#ifdef STARDESTR_EXT
if(m_act[i] == 0.0) dt_new = dt_max;
#endif
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;
// 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<n_act; i++)
{
if( dt_act[i] < min_dt ) min_dt = dt_act[i];
} /* i */
dt_act[0] = min_dt;
dt_act[1] = min_dt;
#else
#ifdef ADD_BH1
min_dt = dt_act[0];
for(i=1; i<n_act; i++)
{
if( dt_act[i] < min_dt ) min_dt = dt_act[i];
} /* i */
dt_act[0] = min_dt;
#endif // ADD_BH1
#endif // ADD_BH2
#ifdef GMC
min_dt = dt_act[0];
for(i=1; i<n_act; i++)
{
if( dt_act[i] < min_dt ) min_dt = dt_act[i];
} /* i */
dt_gmc = min_dt;
#endif // GMC
#ifdef GMC2111
min_dt = dt_act[0];
for(i=1; i<n_act; i++)
{
if( dt_act[i] < min_dt ) min_dt = dt_act[i];
} /* i */
for(i=1; i<n_act; i++)
{
if(m_act[i] > 0.1) dt_act[i] = min_dt;
} /* i */
#endif // GMC2
#ifdef GMC2222
/* define the max. dt for the GMC... */
for(i=1; i<n_act; i++)
{
if(m_act[i] > 0.1) dt_act[i] = dt_max;
} /* i */
#endif // GMC2
#ifdef BBH_INF
if(myRank == rootRank)
{
out = fopen("bbh.inf","a");
i_bh1 = 0;
i_bh2 = 1;
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<n_act; i++)
{
tmp_r2 = SQR(x_act[i][0] - x_bbhc[0]) + SQR(x_act[i][1] - x_bbhc[1]) + SQR(x_act[i][2] - x_bbhc[2]);
iii = ind_act[i];
// if( (tmp_r2 < DR2*R_INF2) )
if( (tmp_r2 < SEMI_a2*R_INF2) )
{
if( (inf_event[iii] == 0) )
{
// printf("INF1: %.6E %.16E %07d %07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .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, i, ind_act[i],
// 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],
// sqrt(tmp_r2),
/// m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
// dt_act[i]);
// fflush(stdout);
fprintf(out,"INF1 %.6E %.16E %07d %07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .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, i, ind_act[i],
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],
sqrt(tmp_r2),
m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
dt_act[i]);
inf_event[iii] = 1;
}
}
else
{
if( (inf_event[iii] == 1) )
{
// printf("INF2: %.6E %.16E %07d %07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .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, i, ind_act[i],
// 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],
// sqrt(tmp_r2),
// m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
// dt_act[i]);
// fflush(stdout);
fprintf(out,"INF2 %.6E %.16E %07d %07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .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, i, ind_act[i],
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],
sqrt(tmp_r2),
m_act[i], x_act[i][0], x_act[i][1], x_act[i][2], v_act[i][0], v_act[i][1], v_act[i][2], pot_act[i],
dt_act[i]);
}
inf_event[iii] = 0;
} /* if(tmp_r2 < DR2*R_INF2) */
} /* i */
fclose(out);
} /* if(myRank == rootRank) */
#endif
/* Return back the new coordinates + etc... of active part. to the global data... */
for(i=0; i<n_act; i++)
{
iii = ind_act[i];
m[iii] = m_act[i];
for(k=0;k<3;k++)
{
x[iii][k] = x_act[i][k];
v[iii][k] = v_act[i][k];
} /* k */
t[iii] = t_act[i];
dt[iii] = dt_act[i];
pot[iii] = pot_act[i];
#ifdef EXTPOT
pot_ext[iii] = pot_act_ext[i];
#endif
for(k=0;k<3;k++)
{
a[iii][k] = a_act[i][k];
adot[iii][k] = adot_act[i][k];
} /* k */
/*
if(myRank == rootRank)
{
tmp_r = sqrt( SQR(x[iii][0]) + SQR(x[iii][1]) + SQR(x[iii][2]) );
if(tmp_r < 1e-8)
{
printf("SOS: TS = %.8E \t t = %.8E \t i = %06d \t R = %.8E \n", Timesteps, time_cur, iii, tmp_r);
fflush(stdout);
exit(-1);
}
}
*/
} /* i */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_ACT_CORR += (CPU_tmp_user - CPU_tmp_user0);
#endif
/* STARDESTR routine on all the nodes */
#ifdef STARDESTR
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
#ifdef ADD_BH2
i_bh = 0;
star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh);
m_act[i_bh] = m_bh1;
i_bh = 1;
star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
N, ind, m, x, v, pot, a, adot, t, &m_bh2, &num_bh2, i_bh);
m_act[i_bh] = m_bh2;
#else
#ifdef ADD_BH1
i_bh = 0;
star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh);
m_act[i_bh] = m_bh1;
#endif // ADD_BH1
#endif // ADD_BH2
#ifdef LIMITS
for(i=0; i<n_act; i++)
{
if(m_act[i] == 0.0)
{
x_act[i][0] = R_LIMITS + 1.0*my_rand2();
x_act[i][1] = R_LIMITS + 1.0*my_rand2();
x_act[i][2] = R_LIMITS + 1.0*my_rand2();
v_act[i][0] = 1.0E-06*my_rand2();
v_act[i][1] = 1.0E-06*my_rand2();
v_act[i][2] = 1.0E-06*my_rand2();
}
}
#endif // LIMITS
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_STARDESTR += (CPU_tmp_user - CPU_tmp_user0);
#endif
#endif // STARDESTR
#ifdef STARDESTR_EXT
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
star_destr_ext(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act,
N, m, x, v, a, adot, t, &m_bh, &num_bh);
// star_destr_ext();
#ifdef LIMITS
for(i=0; i<n_act; i++)
{
if(m_act[i] == 0.0)
{
x_act[i][0] = R_LIMITS + 1.0*my_rand2();
x_act[i][1] = R_LIMITS + 1.0*my_rand2();
x_act[i][2] = R_LIMITS + 1.0*my_rand2();
v_act[i][0] = 1.0E-06*my_rand2();
v_act[i][1] = 1.0E-06*my_rand2();
v_act[i][2] = 1.0E-06*my_rand2();
}
} /* i */
#endif // LIMITS
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_STARDESTR += (CPU_tmp_user - CPU_tmp_user0);
#endif
#endif // STARDESTR_EXT
#if defined(STARDESTR) || defined(STARDESTR_EXT)
// printf("%.0f \t %08d \n", Timesteps+1.0, n_act);
// fflush(stdout);
/* Return back the new coordinates + etc... of active part. to the global data... */
for(i=0; i<n_act; i++)
{
iii = ind_act[i];
m[iii] = m_act[i];
for(k=0;k<3;k++)
{
x[iii][k] = x_act[i][k];
v[iii][k] = v_act[i][k];
} /* k */
t[iii] = t_act[i];
dt[iii] = dt_act[i];
pot[iii] = pot_act[i];
#ifdef EXTPOT
pot_ext[iii] = pot_act_ext[i];
#endif
for(k=0;k<3;k++)
{
a[iii][k] = a_act[i][k];
adot[iii][k] = adot_act[i][k];
} /* k */
} /* i */
#endif // #if defined(STARDESTR) || defined(STARDESTR_EXT)
/* load the new values for active particles to the local GRAPE's */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
for(j=0; j<n_act; j++)
{
#ifdef ETICS_CEP
if (ind_act[j] == grapite_cep_index) grapite_update_cep(t_act[j], x_act[j], v_act[j], a_act[j], adot_act[j]); // All ranks should do it.
#endif
cur_rank = ind_act[j]/n_loc;
if(myRank == cur_rank)
{
jjj = ind_act[j] - myRank*n_loc;
for(k=0;k<3;k++)
{
a2by18[k] = 0.0; a1by6[k] = adot_act[j][k]*over6; aby2[k] = a_act[j][k]*over2;
} /* k */
#ifdef SAPPORO
g6_set_j_particle_(&clusterid, &jjj, &ind_act[j], &t_act[j], &dt_act[j], &m_act[j], a2by18, a1by6, aby2, v_act[j], x_act[j]);
#else
g6_set_j_particle(clusterid, jjj, ind_act[j], t_act[j], dt_act[j], m_act[j], a2by18, a1by6, aby2, v_act[j], x_act[j]);
#endif
} /* if(myRank == cur_rank) */
} /* j */
#ifdef AMD
g6_flush_jp_buffer(clusterid);
#endif
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_ACT_LOAD += (CPU_tmp_user - CPU_tmp_user0);
#endif
/* Current time set to min_t */
time_cur = min_t;
Timesteps += 1.0;
n_act_sum += n_act;
n_act_distr[n_act-1]++;
/*
for(i=1; i<N+1; i++)
{
if(i == n_act) n_act_distr[i-1]++;
}
*/
/* STEVOL routine on all the nodes */
#ifdef STEVOL
#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<N; i++)
{
#ifdef GMC2
if(event[i] == 0) m[i] = m_curr(m0[i], ZZZ[i], t0[i], time_cur);
// if( m[i] < 0.1 ) m[i] = m_curr(m0[i], ZZZ[i], t0[i], time_cur); // exclude the GMC's
#else
if(event[i] == 0) m[i] = m_curr(m0[i], ZZZ[i], t0[i], time_cur);
// m[i] = m_curr(m0[i], ZZZ[i], t0[i], time_cur);
#endif
} /* i */
/*
mcm = 0.0;
for(i=0; i<N; i++) mcm += m[i];
printf("%.4E % .6E \n", time_cur, mcm);
fflush(stdout);
*/
t_stevol += dt_stevol;
} /* if(time_cur >= 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<N; i++)
{
if( m[i] != 0.0 ) m[i] = m_curr(i, m0[i], ZZZ[i], t0[i], time_cur);
} /* i */
/*
mcm = 0.0;
for(i=0; i<N; i++) mcm += m[i];
printf("%.4E % .6E \n", time_cur, mcm);
fflush(stdout);
*/
t_stevol += dt_stevol;
} /* if(time_cur >= 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<N; j++)
{
cur_rank = ind[j]/n_loc;
if(myRank == cur_rank)
{
jjj = ind[j] - myRank*n_loc;
for(k=0;k<3;k++)
{
a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;
} /* k */
#ifdef SAPPORO
g6_set_j_particle_(&clusterid, &jjj, &ind[j], &t[j], &dt[j], &m[j], a2by18, a1by6, aby2, v[j], x[j]);
#else
g6_set_j_particle(clusterid, jjj, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
#endif
} /* if(myRank == cur_rank) */
} /* j */
#ifdef AMD
g6_flush_jp_buffer(clusterid);
#endif
#endif
if(time_cur >= 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<N; i++)
{
for(k=0;k<3;k++)
{
x[i][k] -= xcm[k];
// v[i][k] -= vcm[k];
} /* k */
} /* i */
#endif
/* write cont data */
write_cont_data();
#ifdef GMC
write_cont_GMC();
#endif
/* 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);
// Timesteps, n_act_sum, 57.0*N*n_act_sum/(DT_TOT - DT_ACT_REDUCE)/1.0E+09);
fclose(out);
#endif
#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<N+1;i++)
{
fprintf(tmp_file,"%06d \t %.6E \t %.2f \n", i, n_act_distr[i-1], 100.0*n_act_distr[i-1]/Timesteps);
}
fclose(tmp_file);
#endif
} /* if(myRank == rootRank) */
/* possible coordinate & velocity limits for ALL particles !!! */
#ifdef LIMITS_ALL
for(i=0; i<N; i++)
{
tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
// tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
tmp_rv = x[i][0]*v[i][0] + x[i][1]*v[i][1] + x[i][2]*v[i][2];
// if( (tmp_r > 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<N; j++)
{
cur_rank = ind[j]/n_loc;
if(myRank == cur_rank)
{
jjj = ind[j] - myRank*n_loc;
for(k=0;k<3;k++)
{
a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;
}
#ifdef SAPPORO
g6_set_j_particle_(&clusterid, &jjj, &ind[j], &t[j], &dt[j], &m[j], a2by18, a1by6, aby2, v[j], x[j]);
#else
g6_set_j_particle(clusterid, jjj, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
#endif
}
}
#ifdef AMD
g6_flush_jp_buffer(clusterid);
#endif
#endif // LIMITS_ALL
#ifdef CMCORR111
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* Broadcast the values of all particles to all processors... */
MPI_Bcast(x, 3*N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(v, 3*N, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
/* Scatter the "local" vectors from "global" */
MPI_Scatter(x, 3*n_loc, MPI_DOUBLE, x_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Scatter(v, 3*n_loc, MPI_DOUBLE, v_loc, 3*n_loc, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
#endif
#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();
// 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<N; i++)
{
tmp_r = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
if( tmp_r > 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<N; j++)
{
cur_rank = ind[j]/n_loc;
if(myRank == cur_rank)
{
jjj = ind[j] - myRank*n_loc;
for(k=0;k<3;k++)
{
a2by18[k] = 0.0; a1by6[k] = adot[j][k]*over6; aby2[k] = a[j][k]*over2;
} /* k */
#ifdef SAPPORO
g6_set_j_particle_(&clusterid, &jjj, &ind[j], &t[j], &dt[j], &m[j], a2by18, a1by6, aby2, v[j], x[j]);
#else
g6_set_j_particle(clusterid, jjj, ind[j], t[j], dt[j], m[j], a2by18, a1by6, aby2, v[j], x[j]);
#endif
} /* if(myRank == cur_rank) */
} /* j */
#ifdef AMD
g6_flush_jp_buffer(clusterid);
#endif
#endif // LIMITS_NEW
t_disk += dt_disk;
} /* if(time_cur >= 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<N+1;i++)
{
fprintf(tmp_file,"%06d \t %.6E \t %.2f \n", i, n_act_distr[i-1], 100.0*n_act_distr[i-1]/Timesteps);
}
fclose(tmp_file);
#endif
} /* if(myRank == rootRank) */
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* Finalize the MPI work */
MPI_Finalize();
return(0);
}
/****************************************************************************/
/****************************************************************************/
/****************************************************************************/
Reference in a new issue View git blame Copy permalink