yohai/phigrape
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity
phigrape/phigrape.cpp
Yohai Meiron 5cb4282be4 Completed HDF5 support
2020-03-31 19:04:18 -04:00

2020 lines
68 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
*****************************************************************************/
#include "double3.h"
#include "config.h"
#include "io.h"
Config *config;
#define NORM // Physical normalization
//#define EXTPOT_BH // BH - usually NB units
//#define EXTPOT_GAL_DEH // Dehnen Galactic - usually physical units
#ifdef ETICS
#include "grapite.h"
// why do we need CEP as a compilaion flag... just have it always on when ETICS is on. IF there is no CEP, there should be a graceful skipping of those operations.
//#define ETICS_CEP
#ifndef ETICS_DTSCF
#error "ETICS_DTSCF must be defined"
#endif
#endif
#define TIMING
#define ETA_S_CORR 4.0
#define ETA_BH_CORR 4.0
#define DTMAXPOWER -3.0
#define DTMINPOWER -36.0
/*
-3.0 0.125
-4.0 0.0625
-5.0 0.03125
-7.0 ~1e-2
-10.0 ~1e-3
.............
-20.0 ~1e-6
-23.0 ~1e-7
-26.0 ~1e-8
-30.0 ~1e-9
*/
//#define ACT_DEF_LL
#if defined(ACT_DEF_LL) && defined(ACT_DEF_GRAPITE)
#error "Contradicting preprocessor flags!"
#endif
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <algorithm>
#include <numeric>
#include <stdexcept>
#define G6_NPIPE 2048
#include "grape6.h"
#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)
double L[3]; // needed in pn_bh_spin.c
// Needed for things related to BHs
#include "debug.h"
// int ind_sort[N_MAX];
// double var_sort[N_MAX];
double CPU_time_real0, CPU_time_user0, CPU_time_syst0;
double CPU_time_real, CPU_time_user, CPU_time_syst;
#ifdef TIMING
double CPU_tmp_real0, CPU_tmp_user0, CPU_tmp_syst0;
double CPU_tmp_real, CPU_tmp_user, CPU_tmp_syst;
double DT_TOT,
DT_ACT_DEF1, DT_ACT_DEF2, DT_ACT_DEF3, DT_ACT_PRED,
DT_ACT_GRAV, DT_EXT_GRAV,
DT_GMC_GRAV, DT_GMC_GMC_GRAV, DT_EXT_GMC_GRAV,
DT_ACT_CORR, DT_ACT_LOAD,
DT_STEVOL, DT_STARDISK, DT_STARDESTR;
double DT_ACT_REDUCE;
#endif
/* some local settings for G6a board's */
int clusterid, ii, nn, numGPU;
int npipe=G6_NPIPE, index_i[G6_NPIPE];
double h2_i[G6_NPIPE], p_i[G6_NPIPE];
double3 x_i[G6_NPIPE], v_i[G6_NPIPE],
a_i[G6_NPIPE], jerk_i[G6_NPIPE];
int new_tunit=51, new_xunit=51;
double ti=0.0;
double3 a2by18, a1by6, aby2;
/* normalization... */
double m_norm, r_norm, v_norm, t_norm;
double eps_BH=0.0;
/* external potential... */
#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_BH
double r2, rv_ij,
m_bh, b_bh, eps_bh; // NOTE different from eps_BH
// NOTE there are other m_bh and b_bh defined outside this ifdef block
// NOTE this is a mess. Looks like eps_bh is never used, and the m_bh from outside the block is never used.
#endif
// double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3];
double3 x_bh1, x_bh2, v_bh1, v_bh2;
// double pot_bh1, a_bh1[3], adot_bh1[3],
// pot_bh2, a_bh2[3], adot_bh2[3];
double pot_bh1, pot_bh2;
double3 a_bh1, adot_bh1, a_bh2, adot_bh2;
//double eps_BH = 0.0;
#include "n_bh.c"
double3 a_pn1[7], adot_pn1[7], a_pn2[7], adot_pn2[7];
#include "pn_bh_spin.c"
#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;
gettimeofday(&tv, NULL);
*time_real = tv.tv_sec + 1.0E-06 * tv.tv_usec;
*time_user = *time_real;
}
void write_snap_data(char out_fname[], int diskstep, int N, double time_cur, int ind[], double m[], double3 x[], double3 v[])
{
auto out = fopen(out_fname, "w");
fprintf(out,"%06d \n", diskstep);
fprintf(out,"%07d \n", N);
fprintf(out,"%.10E \n", time_cur);
for (int 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]);
}
fclose(out);
}
void write_bh_data(double time_cur, double m[], double3 x[], double3 v[], double pot[], double3 a[], double3 adot[], double dt[])
{
if (config->live_smbh_count == 2) {
auto out = fopen("bh.dat", "a");
for (int i=0; i < 2; i++) {
double3 *a_pn, *adot_pn, a_bh, adot_bh;
double pot_bh;
if (i==0) {
a_pn = a_pn1;
adot_pn = adot_pn1;
pot_bh = pot_bh1;
a_bh = a_bh1;
adot_bh = adot_bh1;
} else {
a_pn = a_pn2;
adot_pn = adot_pn2;
pot_bh = pot_bh2;
a_bh = a_bh2;
adot_bh = adot_bh2;
}
if (config->live_smbh_custom_eps >= 0) {
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 ",
time_cur, m[i],
x[i][0], x[i][1], x[i][2], x[i].norm(),
v[i][0], v[i][1], v[i][2], v[i].norm(),
pot[i],
a[i][0], a[i][1], a[i][2], a[i].norm(),
adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
dt[i],
pot_bh,
a_bh[0], a_bh[1], a_bh[2], a_bh.norm(),
adot_bh[0], adot_bh[1], adot_bh[2], adot_bh.norm());
if (config->binary_smbh_pn) {
fprintf(out, "\t");
for (int pn_idx=0; pn_idx < 7; pn_idx++) {
fprintf(out, "\t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E ", a_pn[pn_idx][0], a_pn[pn_idx][1], a_pn[pn_idx][2], a_pn[pn_idx].norm(), adot_pn[pn_idx][0], adot_pn[pn_idx][1], adot_pn[pn_idx][2], adot_pn[pn_idx].norm());
}
}
fprintf(out, "\n");
} 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], x[i].norm(),
v[i][0], v[i][1], v[i][2], v[i].norm(),
pot[i],
a[i][0], a[i][1], a[i][2], a[i].norm(),
adot[i][0], adot[i][1], adot[i][2], adot[i].norm(),
dt[i]);
}
}
fprintf(out, "\n");
fclose(out);
} else if (config->live_smbh_count == 1) {
auto out = fopen("bh.dat", "a");
double tmp_r = sqrt( SQR(x[0][0]) + SQR(x[0][1]) + SQR(x[0][2]) );
double tmp_v = sqrt( SQR(v[0][0]) + SQR(v[0][1]) + SQR(v[0][2]) );
double tmp_a = sqrt( SQR(a[0][0]) + SQR(a[0][1]) + SQR(a[0][2]) );
double tmp_adot = sqrt( SQR(adot[0][0]) + SQR(adot[0][1]) + SQR(adot[0][2]) );
fprintf(out,"%.16E \t %.8E \t % .16E % .16E % .16E \t %.16E \t % .16E % .16E % .16E \t %.16E \t % .8E \t % .8E % .8E % .8E \t %.8E \t % .8E % .8E % .8E \t %.8E \t %.8E \n",
time_cur, m[0],
x[0][0], x[0][1], x[0][2], tmp_r,
v[0][0], v[0][1], v[0][2], tmp_v,
pot[0],
a[0][0], a[0][1], a[0][2], tmp_a,
adot[0][0], adot[0][1], adot[0][2], tmp_adot,
dt[0]);
fprintf(out,"\n");
fclose(out);
}
}
void write_bh_nb_data(double time_cur, int N, double m[], double3 x[], double3 v[])
{
int nb = config->live_smbh_neighbor_number;
int ind_sort[N_MAX];
double var_sort[N_MAX];
auto out = fopen("bh_neighbors.dat", "a");
/* 1st BH */
for (int i_bh=0; i_bh < config->live_smbh_count; i_bh++) {
for (int i=0; i<N; i++) var_sort[i] = (x[i]-x[i_bh]).norm();
std::iota(ind_sort, ind_sort+N, 0);
std::partial_sort(ind_sort, ind_sort + nb, ind_sort + N, [&](int i, int j) {return var_sort[i] < var_sort[j];});
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 (int j=1; j < nb; j++) {
int i = ind_sort[j];
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], (x[i]-x[i_bh]).norm(),
v[i][0], v[i][1], v[i][2], (v[i]-v[i_bh]).norm());
}
fprintf(out,"\n");
}
fprintf(out, "\n"); // this is redundant
fclose(out);
}
class External_gravity {
public:
void apply(const int n_act, const double3 x[], const double3 v[], double pot[], double3 a[], double3 adot[])
{
for (int i=0; i<n_act; i++) {
this->set_coordinates(x[i], v[i]);
this->calc_gravity();
pot[i] += potential;
a[i] += acceleration;
adot[i] += jerk;
}
}
virtual void calc_gravity() = 0;
bool is_active;
protected:
double potential;
double3 acceleration, jerk;
double3 x, v;
void set_coordinates(double3 x, double3 v)
{
this->x = x;
this->v = v;
}
};
class Plummer : public External_gravity {
public:
Plummer(double m, double b) : m(m), b(b) {is_active=(m>0);}
private:
void calc_gravity()
{
double r2 = SQR(b);
r2 += x.norm2();
double r = sqrt(r2);
double rv_ij = v*x;
double tmp = m / r;
potential = -tmp;
tmp /= r2;
acceleration = - tmp * x;
jerk = - tmp * (v - 3*rv_ij * x/r2);
}
double m, b;
};
class Miyamoto_Nagai : public External_gravity {
public:
Miyamoto_Nagai(double m, double a, double b) : m(m), a(a), b(b) {is_active=(m>0);}
private:
void calc_gravity()
{
double x_ij=x[0], y_ij=x[1], z_ij=x[2];
double vx_ij=v[0], vy_ij=v[1], vz_ij=v[2];
auto z2_tmp = z_ij*z_ij + SQR(b);
auto z_tmp = sqrt(z2_tmp);
auto r2_tmp = x_ij*x_ij + y_ij*y_ij + SQR(z_tmp + a);
auto r_tmp = sqrt(r2_tmp);
potential = - m / r_tmp;
auto tmp = m / (r2_tmp*r_tmp);
acceleration[0] = - tmp * x_ij;
acceleration[1] = - tmp * y_ij;
acceleration[2] = - tmp * z_ij * (z_tmp + a)/z_tmp;
tmp = m / (z_tmp*r2_tmp*r2_tmp*r_tmp);
jerk[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));
jerk[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));
jerk[2] = tmp * (- vz_ij*(z_tmp + a)*(x_ij*x_ij*(z2_tmp*z_tmp + a*SQR(b)) +
y_ij*y_ij*(z2_tmp*z_tmp + a*SQR(b)) -
(2.0*a*(SQR(z_ij) - SQR(b))*z_tmp +
2.0*SQR(z_ij*z_ij) +
SQR(b)*SQR(z_ij) -
SQR(b)*(SQR(a) + SQR(b))))
+ 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a)
+ 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a)) / z2_tmp;
}
double m, a, b;
};
class Logarithmic_halo : public External_gravity {
public:
Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0);}
private:
void calc_gravity()
{
auto r2 = x.norm2();
auto rv_ij = 2.0*(v*x);
auto r2_r2_halo = (r2 + r2_halo);
potential = - 0.5*v2_halo * log(1.0 + r2/r2_halo);
auto tmp = v2_halo/r2_r2_halo;
acceleration = - tmp * x;
tmp /= (r2 + r2_halo);
jerk = tmp * (rv_ij * x - r2_r2_halo * v);
}
double v2_halo, r2_halo;
};
void calc_self_grav(double t, double eps2, double &g6_calls, int n_loc,
int n_act, int ind_act[],
double3 x_act_new[], double3 v_act_new[],
double pot_act_tmp[],
double3 a_act_tmp[],
double3 adot_act_tmp[],
double h2_i[])
{
/* calc the new grav for the active particles */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
g6_set_ti(clusterid, t);
int ni = n_act; // TODO why is this needed?
/* define the local phi, a, adot for these active particles */
for (int i=0; i<ni; i+=npipe) {
nn = npipe;
if (ni-i < npipe) nn = ni - i;
for (ii=0; ii<nn; ii++) {
h2_i[ii] = eps2; // TODO This should be a global or something
} /* ii */
//TODO any way we can clean up this ugly casting?
g6calc_firsthalf(clusterid, n_loc, nn, ind_act+i, (double(*)[3])x_act_new+i, (double(*)[3])v_act_new+i, (double(*)[3])a_act_tmp+i, (double(*)[3])adot_act_tmp+i, pot_act_tmp+i, eps2, h2_i);
g6calc_lasthalf( clusterid, n_loc, nn, ind_act+i, (double(*)[3])x_act_new+i, (double(*)[3])v_act_new+i, eps2, h2_i, (double(*)[3])a_act_tmp+i, (double(*)[3])adot_act_tmp+i, pot_act_tmp+i);
g6_calls++;
} /* 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(std::vector<External_gravity*> &external_gravity_components, int n_act, double3 *x_act_new, double3 *v_act_new, double *pot_act_ext, double3 *a_act_new, double3* adot_act_new)
{
/* Define the external potential for all active particles on all nodes */
int ni = n_act; // TODO redundant?
std::fill(pot_act_ext, pot_act_ext+n_act, 0.);
for (auto component : external_gravity_components) {
if (component->is_active)
component->apply(n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new);
}
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#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 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 = 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 TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
#endif
}
void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, double& rcm_sum, double& vcm_sum, double& E_tot_0, double& E_tot_corr_0, double& E_tot_corr_sd_0, int &skip_con, int N, double m[], double3 x[], double3 v[], double pot[], double pot_ext[])
{
// TODO maybe use static variables here for the previous step energy?
double E_pot = 0.0;
for (int i=0; i<N; i++) E_pot += m[i]*pot[i];
E_pot *= 0.5;
double E_kin = 0.0;
for (int 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;
double E_pot_ext = 0.0;
for (int i=0; i<N; i++) E_pot_ext += m[i]*pot_ext[i];
if (timesteps == 0.0) { //TODO what is this thing?
rcm_sum = 0.0;
vcm_sum = 0.0;
}
double mcm = 0.0;
double rcm_mod = 0.0;
double vcm_mod = 0.0;
double3 xcm = {0, 0, 0};
double3 vcm = {0, 0, 0};
for (int i=0; i<N; i++) {
mcm += m[i];
xcm += m[i] * x[i];
vcm += m[i] * v[i];
} /* i */
xcm /= mcm;
vcm /= mcm;
rcm_mod = xcm.norm();
vcm_mod = vcm.norm();
rcm_sum += rcm_mod;
vcm_sum += vcm_mod;
double3 mom = {0, 0, 0};
for (int 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);
double E_tot = E_pot + E_kin + E_pot_ext;
double E_tot_corr = E_tot /* + E_corr */; // TODO what is this?
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;
skip_con = 1;
} else {
/* read the E_tot_0 + etc... */
if (skip_con == 0) {
auto inp = fopen("contr.con","r");
double tmp;
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;
}
}
double DE_tot = E_tot - E_tot_0;
/* This is the only output to screen */
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);
fflush(stdout);
auto 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;
E_tot_corr_0 = E_tot_corr;
E_tot_corr_sd_0 = 0;
}
int main(int argc, char *argv[])
{
int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank,
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,
eta_s, eta, eta_bh,
E_tot_0, E_tot_corr_0, E_tot_corr_sd_0,
rcm_sum=0.0, vcm_sum=0.0,
eps=0.0, eps2,
a2_mod, adot2_mod,
dt_tmp, dt2half, dt3over6, dt4over24, dt5over120,
dtinv, dt2inv, dt3inv,
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;
double3 xcm, vcm, mom,
xdc, vdc,
a2, a3, a2dot1;
char processor_name[MPI_MAX_PROCESSOR_NAME],
inp_fname[30];
/* global variables */
int N, ind[N_MAX];
double m[N_MAX], pot[N_MAX], t[N_MAX], dt[N_MAX];
double3 x[N_MAX], v[N_MAX], a[N_MAX], adot[N_MAX];
/* local variables */
int n_loc, ind_loc[N_MAX_loc];
double m_loc[N_MAX_loc], pot_loc[N_MAX_loc], t_loc[N_MAX_loc], dt_loc[N_MAX_loc];
double3 x_loc[N_MAX_loc], v_loc[N_MAX_loc],
a_loc[N_MAX_loc], adot_loc[N_MAX_loc];
/* data for active particles */
int n_act, ind_act[N_MAX];
double m_act[N_MAX],
pot_act[N_MAX], t_act[N_MAX], dt_act[N_MAX],
pot_act_new[N_MAX],
pot_act_tmp[N_MAX];
double3 x_act[N_MAX], v_act[N_MAX],
a_act[N_MAX], adot_act[N_MAX],
x_act_new[N_MAX], v_act_new[N_MAX],
a_act_new[N_MAX], adot_act_new[N_MAX],
a_act_tmp[N_MAX], adot_act_tmp[N_MAX];;
FILE *out;
double C_NB;
int usedOrNot[7];
double pot_ext[N_MAX], pot_act_ext[N_MAX]; // for all EXTPOT
int i_bh1, i_bh2;
double m_bh1, m_bh2;
int inf_event[N_MAX];
double DR2, tmp_r2;
double3 x_bbhc, v_bbhc;
double DV2, EB, SEMI_a, SEMI_a2;
double s_bh1[3] = {0.0, 0.0, 1.0};
double s_bh2[3] = {0.0, 0.0, 1.0};
/* INIT the rand() !!! */
srand(19640916); /* it is just my birthday :-) */
/* 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);
/* read the input parameters to the rootRank */
config = new Config("phigrape.conf");
C_NB = config->pn_c;
std::copy(config->pn_usage.begin(), config->pn_usage.end(), usedOrNot);
if (is_hdf5(config->input_file_name)) {
#ifndef HAS_HDF5
fprintf(stderr, "ERROR: input file is in HDF5 format, but the code was compiled without HDF5 support\n")
return -1;
#endif
h5_read(config->input_file_name, &diskstep, &N, &time_cur, m, x, v);
}
else
ascii_read(config->input_file_name, &diskstep, &N, &time_cur, m, x, v);
std::iota(ind, ind+N, 0);
if (myRank == rootRank) {
//TODO move it out of (myRank == rootRank) so you don't need to communicate them.
eps = config->eps;
t_end = config->t_end;
dt_disk = config->dt_disk;
dt_contr = config->dt_contr;
dt_bh = config->dt_bh;
eta = config->eta;
strcpy(inp_fname, config->input_file_name.c_str());
if (config->binary_smbh_influence_sphere_output) for (int i=0; i<N; i++) inf_event[i] = 0; // WARNING N wasn't set yet!
/*
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_neighbors.dat)
eta : parameter for timestep determination
inp_data : name of the input file (data.inp)
*/
/* Read the data for EXTernal POTential */
#ifdef EXTPOT_BH
fscanf(inp,"%lE %lE", &m_bh, &b_bh);
#endif
#ifdef EXTPOT_GAL_DEH
fscanf(inp,"%lE %lE %lE", &m_ext, &r_ext, &g_ext);
#endif
printf("\n");
printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc);
printf("\n");
printf("N = %07d \t eps = %.6E \n", N, eps);
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
printf("External Potential: \n");
printf("\n");
#ifdef EXTPOT_BH
printf("m_bh = %.4E b_bh = %.4E \n", m_bh, b_bh);
#endif
#ifdef EXTPOT_GAL_DEH
printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n", m_ext, r_ext, g_ext);
#endif
if ((diskstep == 0) && (time_cur == 0.0)) {
out = fopen("contr.dat", "w");
fclose(out);
#ifdef TIMING
out = fopen("timing.dat", "w");
fclose(out);
#endif
if (config->live_smbh_output && (config->live_smbh_count > 0)) {
out = fopen("bh.dat", "w");
fclose(out);
}
if ((config->live_smbh_neighbor_output) && (config->live_smbh_count > 0)) {
out = fopen("bh_neighbors.dat", "w");
fclose(out);
}
if (config->binary_smbh_influence_sphere_output) {
out = fopen("bbh_inf.dat", "w");
fclose(out);
}
} else { // if (diskstep == 0)
} // if (diskstep == 0)
get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0);
} /* if (myRank == rootRank) */
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* Broadcast all useful values to all processors... */
MPI_Bcast(&N, 1, MPI_INT, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&eps, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&eta, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&t_end, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&dt_disk, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&dt_contr, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&dt_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&time_cur, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
#ifdef NORM
MPI_Bcast(&m_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&r_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&v_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&t_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
#endif // NORM
#ifdef EXTPOT_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
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
double normalization_mass=1, normalization_length=1, normalization_velocity=1;
if (config->ext_units_physical) {
normalization_mass = 1/config->unit_mass;
normalization_length = 1000/config->unit_length;
normalization_velocity = 1.52484071426404437233e+01*sqrt(config->unit_length/config->unit_mass);
}
Plummer ext_bulge(config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length);
Miyamoto_Nagai ext_disk(config->ext_m_disk*normalization_mass, config->ext_a_disk*normalization_length, config->ext_b_disk*normalization_length);
Plummer ext_halo_plummer(config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length);
Logarithmic_halo ext_log_halo(config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length);
std::vector<External_gravity*> external_gravity_components;
external_gravity_components.push_back(&ext_bulge);
external_gravity_components.push_back(&ext_disk);
external_gravity_components.push_back(&ext_halo_plummer);
external_gravity_components.push_back(&ext_log_halo);
if (ext_bulge.is_active) printf("m_bulge = %.4E b_bulge = %.4E\n", config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length);
if (ext_disk.is_active) printf("m_disk = %.4E a_disk = %.4E b_disk = %.4E\n", config->ext_m_disk*normalization_mass, config->ext_a_disk*normalization_length, config->ext_b_disk*normalization_length);
if (ext_halo_plummer.is_active) printf("m_halo = %.4E b_halo = %.4E\n", config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length);
if (ext_log_halo.is_active) printf("v_halo = %.6E r_halo = %.6E \n", config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length);
printf("\n");
fflush(stdout);
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) */
std::fill(t, t+N, time_cur);
std::fill(dt, dt+N, 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);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* 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 MPI_OVERRIDE
numGPU = 1; // TODO get this from config file
clusterid = myRank % numGPU;
#else
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);
#endif
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);
#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 (int j=0; j<nj; j++) {
a2by18 = {0, 0, 0};
a1by6 = {0, 0, 0};
aby2 = {0, 0, 0};
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]);
} /* j */
#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 (int i=0; i<n_act; i++) {
ind_act[i] = ind[i];
iii = ind_act[i];
m_act[i] = m[iii];
x_act[i] = x[iii];
v_act[i] = v[iii];
t_act[i] = t[iii];
dt_act[i] = dt[iii];
} /* i */
// NOTE this is where calc_self_grav_zero() used to be.
calc_self_grav(time_cur, eps2, g6_calls, n_loc,
n_act, ind_act,
x_act, v_act,
pot_act_tmp,
a_act_tmp,
adot_act_tmp,
h2_i);
/* 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);
if (config->live_smbh_count == 2) {
i_bh1 = 0;
i_bh2 = 1;
if (config->live_smbh_custom_eps >= 0) {
m_bh1 = m[i_bh1];
m_bh2 = m[i_bh2];
x_bh1 = x[i_bh1];
v_bh1 = v[i_bh1];
x_bh2 = x[i_bh2];
v_bh2 = v[i_bh2];
// calculate and "minus" the BH <-> BH _softened_ pot, acc & jerk
tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
m_bh2, x_bh2, v_bh2,
eps,
&pot_bh1, a_bh1, adot_bh1,
&pot_bh2, a_bh2, adot_bh2);
pot[i_bh1] -= pot_bh1;
pot[i_bh2] -= pot_bh2;
a[i_bh1] -= a_bh1;
a[i_bh2] -= a_bh2;
adot[i_bh1] -= adot_bh1;
adot[i_bh2] -= adot_bh2;
// calculate and "plus" the new BH <-> BH _unsoftened_ pot, acc, jerk
tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
m_bh2, x_bh2, v_bh2,
config->live_smbh_custom_eps,
&pot_bh1, a_bh1, adot_bh1,
&pot_bh2, a_bh2, adot_bh2);
pot[i_bh1] += pot_bh1;
pot[i_bh2] += pot_bh2;
a[i_bh1] += a_bh1;
a[i_bh2] += a_bh2;
adot[i_bh1] += adot_bh1;
adot[i_bh2] += adot_bh2;
}
if (config->binary_smbh_pn) {
// calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
dt_bh_tmp = dt[0];
tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1,
m_bh2, x_bh2, v_bh2, s_bh2,
C_NB, dt_bh_tmp, usedOrNot,
(double(*)[3])a_pn1, (double(*)[3])adot_pn1,
(double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank);
a[i_bh1] += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6];
a[i_bh2] += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6];
adot[i_bh1] += adot_pn1[1] + adot_pn1[2] + adot_pn1[3] + adot_pn1[4] + adot_pn1[5] + adot_pn1[6];
adot[i_bh2] += adot_pn2[1] + adot_pn2[2] + adot_pn2[3] + adot_pn2[4] + adot_pn2[5] + adot_pn2[6];
if (myRank == rootRank) {
if (tmp_i == 505) {
printf("PN RSDIST: %.8E \t %.8E \n", timesteps, time_cur);
fflush(stdout);
exit(-1);
}
}
}
}
calc_ext_grav(external_gravity_components, n_act, x_act, v_act, pot_ext, a, adot);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* Energy control... */
if (myRank == rootRank) {
energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot, pot_ext);
} /* 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 (int 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;
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[i] = dt_tmp;
if (config->dt_min_warning && (myRank == 0)) {
if (dt[i] == dt_min) {
printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[i], ind[i]);
fflush(stdout);
}
}
} /* i */
if (config->live_smbh_count > 0) {
double min_dt = *std::min_element(dt, dt+N);
for (int i=0; i<config->live_smbh_count; i++) dt[i] = min_dt;
}
/* 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 (int j=0; j<nj; j++) {
a2by18 = {0, 0, 0};
a1by6 = adot_loc[j]*(1./6.);
aby2 = a_loc[j]*0.5;
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]);
} /* j */
if (myRank == rootRank) {
/* Write BH data... */
if (config->live_smbh_output) write_bh_data(time_cur, m, x, v, pot, a, adot, dt);
/* Write BH NB data... */
if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur, N, m, x, v);
} /* if (myRank == rootRank) */
/* 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);
} /* if (myRank == rootRank) */
timesteps = 0.0;
n_act_sum = 0.0;
for (int i=1; i<N+1; i++) { // TODO why like this?
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 */
while (time_cur <= t_end) {
/* Define the minimal time and the active particles on all the nodes (exclude the ZERO masses!!!) */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
#ifdef ACT_DEF_GRAPITE
min_t_loc = grapite_get_minimum_time();
#else
min_t_loc = t[0]+dt[0];
for (int 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;
for (int i=0; i<N; i++) {
if (t[i]+dt[i] == min_t ) {
ind_act[n_act] = i;
n_act++;
}
} /* i */
#endif // ACT_DEF_GRAPITE
// TODO deal with it below
#ifdef ACT_DEF_GRAPITE
#error please fix here
#endif
#if defined(ACT_DEF_GRAPITE) && (defined(ADD_BH1) || defined(ADD_BH2))
#ifdef ADD_BH1
#define ACT_DEF_GRAPITE_NUMBH 1
#else
#define ACT_DEF_GRAPITE_NUMBH 2
#endif
int acf_def_grapite_bh_count = 0;
for (i=0; i<n_act; i++) {
if (ind_act[i]==0) {
i_bh1 = i;
acf_def_grapite_bh_count++;
}
#ifdef ADD_BH2
else if (ind_act[i]==1) {
i_bh2 = i;
acf_def_grapite_bh_count++;
}
#endif
if (acf_def_grapite_bh_count==ACT_DEF_GRAPITE_NUMBH) break;
}
if (i==n_act) {
fprintf(stderr, "ERROR: black holes were not found in the active particle list");
exit(1);
}
#else
if (config->live_smbh_count > 0) {
i_bh1 = 0;
i_bh2 = 1;
}
#endif
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_ACT_DEF2 += (CPU_tmp_user - CPU_tmp_user0);
#endif
#ifdef TIMING
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
for (int i=0; i<n_act; i++) {
iii = ind_act[i];
m_act[i] = m[iii];
x_act[i] = x[iii];
v_act[i] = v[iii];
t_act[i] = t[iii];
dt_act[i] = dt[iii];
pot_act[i] = pot[iii];
pot_act_ext[i] = pot_ext[iii];
a_act[i] = a[iii];
adot_act[i] = adot[iii];
} /* i */
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
DT_ACT_DEF3 += (CPU_tmp_user - CPU_tmp_user0);
#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 (int i=0; i<n_act; i++) {
dt_tmp = min_t - t_act[i];
dt2half = 0.5*dt_tmp*dt_tmp;
dt3over6 = (1./3.)*dt_tmp*dt2half;
x_act_new[i] = x_act[i] + v_act[i]*dt_tmp + a_act[i]*dt2half + adot_act[i]*dt3over6;
v_act_new[i] = v_act[i] + a_act[i]*dt_tmp + adot_act[i]*dt2half;
} /* 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
calc_self_grav(min_t, eps2, g6_calls, n_loc,
n_act, ind_act,
x_act_new, v_act_new,
pot_act_tmp,
a_act_tmp,
adot_act_tmp,
h2_i);
/* 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
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
if (config->live_smbh_count == 2) {
if (config->live_smbh_custom_eps >= 0) {
m_bh1 = m_act[i_bh1];
m_bh2 = m_act[i_bh2];
x_bh1 = x_act_new[i_bh1];
v_bh1 = v_act_new[i_bh1];
x_bh2 = x_act_new[i_bh2];
v_bh2 = v_act_new[i_bh2];
// calculate and "minus" the BH <-> BH softened pot, acc & jerk
tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
m_bh2, x_bh2, v_bh2,
eps,
&pot_bh1, a_bh1, adot_bh1,
&pot_bh2, a_bh2, adot_bh2);
pot_act_new[i_bh1] -= pot_bh1;
pot_act_new[i_bh2] -= pot_bh2;
a_act_new[i_bh1] -= a_bh1;
a_act_new[i_bh2] -= a_bh2;
adot_act_new[i_bh1] -= adot_bh1;
adot_act_new[i_bh2] -= adot_bh2;
// calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk
tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1,
m_bh2, x_bh2, v_bh2,
config->live_smbh_custom_eps,
&pot_bh1, a_bh1, adot_bh1,
&pot_bh2, a_bh2, adot_bh2);
pot_act_new[i_bh1] += pot_bh1;
pot_act_new[i_bh2] += pot_bh2;
a_act_new[i_bh1] += a_bh1;
a_act_new[i_bh2] += a_bh2;
adot_act_new[i_bh1] += adot_bh1;
adot_act_new[i_bh2] += adot_bh2;
}
if (config->binary_smbh_pn) {
// calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
dt_bh_tmp = dt[0];
tmp_i = calc_force_pn_BH(m_bh1, x_bh1, v_bh1, s_bh1,
m_bh2, x_bh2, v_bh2, s_bh2,
C_NB, dt_bh_tmp, usedOrNot,
(double(*)[3])a_pn1, (double(*)[3])adot_pn1,
(double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank);
a_act_new[i_bh1] += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6];
a_act_new[i_bh2] += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6];
adot_act_new[i_bh1] += adot_pn1[1] + adot_pn1[2] + adot_pn1[3] + adot_pn1[4] + adot_pn1[5] + adot_pn1[6];
adot_act_new[i_bh2] += adot_pn2[1] + adot_pn2[2] + adot_pn2[3] + adot_pn2[4] + adot_pn2[5] + adot_pn2[6];
if (myRank == rootRank) {
if (tmp_i == 505) {
printf("PN RSDIST: TS = %.8E \t t = %.8E \n", timesteps, time_cur);
fflush(stdout);
exit(-1);
}
}
}
}
calc_ext_grav(external_gravity_components, n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new);
/* 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 (int 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;
double3 a0mia1 = a_act[i]-a_act_new[i];
double3 ad04plad12 = 4.0*adot_act[i] + 2.0*adot_act_new[i];
double3 ad0plad1 = adot_act[i] + adot_act_new[i];
a2 = -6.0*a0mia1*dt2inv - ad04plad12*dtinv;
a3 = 12.0*a0mia1*dt3inv + 6.0*ad0plad1*dt2inv;
x_act_new[i] += dt4over24*a2 + dt5over120*a3;
v_act_new[i] += dt3over6*a2 + dt4over24*a3;
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]);
a2dot1 = a2 + dt_tmp*a3;
a2dot1abs = a2dot1.norm();
a3dot1abs = a3.norm();
if ((config->live_smbh_count > 0) && (ind_act[i] < config->live_smbh_count))
dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
else
dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
if (dt_new < dt_min) dt_tmp = dt_min;
if ((dt_new < dt_tmp) && (dt_new > dt_min)) {
power = log(dt_new)/log(2.0) - 1;
dt_tmp = pow(2.0, (double)power);
}
if ((dt_new > 2.0*dt_tmp) && (fmod(min_t, 2.0*dt_tmp) == 0.0) && (2.0*dt_tmp <= dt_max)) {
dt_tmp *= 2.0;
}
dt_act[i] = dt_tmp;
t_act[i] = min_t;
pot_act[i] = pot_act_new[i];
x_act[i] = x_act_new[i];
v_act[i] = v_act_new[i];
a_act[i] = a_act_new[i];
adot_act[i] = adot_act_new[i];
if (config->dt_min_warning && (myRank == 0)) {
if (dt_act[i] == dt_min) {
printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]);
fflush(stdout);
}
}
} /* i */
/* define the min. dt over all the act. part. and set it also for the BH... */
if (config->live_smbh_count > 0) {
double min_dt = *std::min_element(dt_act, dt_act+n_act);
if (config->live_smbh_count>=1) dt_act[i_bh1] = min_dt;
if (config->live_smbh_count==2) dt_act[i_bh2] = min_dt;
}
if (config->binary_smbh_influence_sphere_output) {
if (myRank == rootRank) {
out = fopen("bbh_inf.dat", "a");
m_bh1 = m_act[i_bh1];
m_bh2 = m_act[i_bh2];
x_bh1 = x_act[i_bh1];
x_bh2 = x_act[i_bh2];
v_bh1 = v_act[i_bh1];
v_bh2 = v_act[i_bh2];
x_bbhc = (m_bh1*x_bh1 + m_bh2*x_bh2)/(m_bh1 + m_bh2);
v_bbhc = (m_bh1*v_bh1 + m_bh2*v_bh2)/(m_bh1 + m_bh2);
DR2 = SQR(x_bh1[0] - x_bh2[0]) + SQR(x_bh1[1] - x_bh2[1]) + SQR(x_bh1[2] - x_bh2[2]);
DV2 = SQR(v_bh1[0] - v_bh2[0]) + SQR(v_bh1[1] - v_bh2[1]) + SQR(v_bh1[2] - v_bh2[2]);
EB = -(m_bh1 + m_bh2) / sqrt(DR2) + 0.5 * DV2;
SEMI_a = -0.5 * (m_bh1 + m_bh2)/EB;
SEMI_a2 = SQR(SEMI_a);
for (int 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 < SEMI_a2*SQR(config->binary_smbh_influence_radius_factor)) {
if (inf_event[iii] == 0) {
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) {
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) */
}
/* Return back the new coordinates + etc... of active part. to the global data... */
for (int i=0; i<n_act; i++) {
iii = ind_act[i];
m[iii] = m_act[i];
x[iii] = x_act[i];
v[iii] = v_act[i];
t[iii] = t_act[i];
dt[iii] = dt_act[i];
pot[iii] = pot_act[i];
pot_ext[iii] = pot_act_ext[i];
a[iii] = a_act[i];
adot[iii] = adot_act[i];
} /* 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
/* 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 (int j=0; j<n_act; j++) { // TODO would be nicer to use i instead of j here
#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;
a2by18 = {0, 0, 0};
a1by6 = adot_act[j]*(1./6.);
aby2 = a_act[j]*0.5;
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]);
} /* if (myRank == cur_rank) */
} /* j */
#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]++;
if (time_cur >= t_bh) {
if (myRank == rootRank) {
/* Write BH data... */
if (config->live_smbh_output) write_bh_data(time_cur, m, x, v, pot, a, adot, dt);
/* Write BH NB data... */
if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur, N, m, x, v);
} /* if (myRank == rootRank) */
t_bh += dt_bh;
} /* if (time_cur >= t_bh) */
if (time_cur >= t_contr) {
if (myRank == rootRank) {
energy_contr(time_cur, timesteps, n_act_sum, g6_calls, rcm_sum, vcm_sum, E_tot_0, E_tot_corr_0, E_tot_corr_sd_0, skip_con, N, m, x, v, pot, pot_ext);
/* write cont data */
if (config->output_hdf5) h5_write("data.con", diskstep, N, time_cur, m, x, v, pot, a, adot, 0, true);
else ascii_write("data.con", diskstep, N, time_cur, m, x, v, 16);
/* possible OUT for timing !!! */
#ifdef TIMING
out = fopen("timing.dat", "a");
DT_TOT = DT_ACT_DEF1 + DT_ACT_DEF2 + DT_ACT_DEF3 + DT_ACT_PRED +
DT_ACT_GRAV + DT_EXT_GRAV + DT_GMC_GRAV +
DT_GMC_GMC_GRAV + DT_EXT_GMC_GRAV +
DT_ACT_CORR + DT_ACT_LOAD +
DT_STEVOL + DT_STARDISK + DT_STARDESTR +
DT_ACT_REDUCE;
fprintf(out,"%.8E \t %.6E \t %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f \t %.3f \t %.8E %.8E %.8E \t %.8E %.8E %.8E \n",
time_cur, DT_TOT,
100.0*DT_ACT_DEF1/DT_TOT, 100.0*DT_ACT_DEF2/DT_TOT, 100.0*DT_ACT_DEF3/DT_TOT, 100.0*DT_ACT_PRED/DT_TOT,
100.0*DT_ACT_GRAV/DT_TOT, 100.0*DT_EXT_GRAV/DT_TOT, 100.0*DT_GMC_GRAV/DT_TOT,
100.0*DT_GMC_GMC_GRAV/DT_TOT, 100.0*DT_EXT_GMC_GRAV/DT_TOT,
100.0*DT_ACT_CORR/DT_TOT, 100.0*DT_ACT_LOAD/DT_TOT,
100.0*DT_STEVOL/DT_TOT, 100.0*DT_STARDISK/DT_TOT, 100.0*DT_STARDESTR/DT_TOT,
100.0*DT_ACT_REDUCE/DT_TOT,
CPU_time_real-CPU_time_real0, CPU_time_user-CPU_time_user0, CPU_time_syst-CPU_time_syst0,
timesteps, n_act_sum, 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09);
fclose(out);
#endif
} /* if (myRank == rootRank) */
#ifdef ETICS_CEP
// We are /inside/ a control step, so all particles must be synchronized; we can safely calculate their density centre. The acceleration and jerk currently in the memory are for the predicted position of the CEP, by calling grapite_calc_center we "correct" the position and velocity, but not the gravity at that point.
// First calculate the DC
grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc);
// Now copy it to the global particle list
memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
memcpy(v[grapite_cep_index], vdc, 3*sizeof(double));
// Now copy it to the local particle list for tha appropriate rank
if (myRank == grapite_cep_index/n_loc) {
memcpy(x_loc[grapite_cep_index - myRank*n_loc], xdc, 3*sizeof(double));
memcpy(v_loc[grapite_cep_index - myRank*n_loc], vdc, 3*sizeof(double));
}
grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
#endif
t_contr += dt_contr;
} /* if (time_cur >= t_contr) */
if (time_cur >= t_disk) {
if (myRank == rootRank) {
diskstep++;
char out_fname[256];
sprintf(out_fname, "%06d", diskstep);
if (config->output_hdf5) h5_write(std::string(out_fname) + ".h5", diskstep, N, time_cur, m, x, v, pot, a, adot, 0, true);
else ascii_write(std::string(out_fname) + ".dat", diskstep, N, time_cur, m, x, v, 10);
// TODO custom precision
} /* if (myRank == rootRank) */
#ifdef ETICS_DUMP
sprintf(out_fname, "coeffs.%06d.%02d.dat", diskstep, myRank);
grapite_dump(out_fname, 2);
#endif
t_disk += dt_disk;
} /* if (time_cur >= t_disk) */
} /* while (time_cur < t_end) */
/* close the local GRAPEs */
g6_close(clusterid);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
MPI_Reduce(&g6_calls, &g6_calls_sum, 1, MPI_DOUBLE, MPI_SUM, rootRank, MPI_COMM_WORLD);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
if (myRank == rootRank) {
/* Write some output for the timestep annalize... */
printf("\n");
printf("timesteps = %.0f Total sum of integrated part. = %.0f g6_calls on all nodes = %.0f \n", timesteps, n_act_sum, g6_calls);
printf("\n");
printf("Real Speed = %.3f GFlops \n", 57.0*N*n_act_sum/(CPU_time_user-CPU_time_user0)/1.0E+09);
fflush(stdout);
} /* if (myRank == rootRank) */
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* Finalize the MPI work */
MPI_Finalize();
return 0;
}
Reference in a new issue View git blame Copy permalink