yohai/phigrape
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Lots of cleanup and started moving the PN out of the main

Browse source
This commit is contained in:
Yohai Meiron 2020-04-18 21:48:46 -04:00
parent 30ae8631a9
commit 2a50f0fc9a
3 changed files with 278 additions and 219 deletions
Download patch file Download diff file Expand all files Collapse all files

51
double3.h
View file

@ -2,81 +2,80 @@
#include <cmath>
struct double3 {
double data[3];
union {
struct {double x, y, z;};
double data[3]; // For legacy access
};
double3() {}
double3(const double x, const double y, const double z)
{
data[0] = x;
data[1] = y;
data[2] = z;
}
: x(x), y(y), z(z) {}
double& operator[](int i) {return data[i];}
const double operator[](int i) const {return data[i];}
double3& operator=(const double3& a)
{
data[0] = a.data[0];
data[1] = a.data[1];
data[2] = a.data[2];
x = a.x;
y = a.y;
z = a.z;
return *this;
}
double3& operator+=(const double3& a)
{
data[0] += a.data[0];
data[1] += a.data[1];
data[2] += a.data[2];
x += a.x;
y += a.y;
z += a.z;
return *this;
}
double3& operator-=(const double3& a)
{
data[0] -= a.data[0];
data[1] -= a.data[1];
data[2] -= a.data[2];
x -= a.x;
y -= a.y;
z -= a.z;
return *this;
}
double3& operator/=(const double& c)
{
data[0] /= c;
data[1] /= c;
data[2] /= c;
x /= c;
y /= c;
z /= c;
return *this;
}
double norm2() const
{
return data[0]*data[0]+data[1]*data[1]+data[2]*data[2];
return x*x + y*y + z*z;
}
double norm() const
{
return sqrt(data[0]*data[0]+data[1]*data[1]+data[2]*data[2]);
return sqrt(x*x + y*y + z*z);
}
operator double*() {return data;}
};
inline double3 operator*(const double& c, const double3& a)
{
return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c);
return double3(a.x*c, a.y*c, a.z*c);
}
inline double3 operator*(const double3& a, const double& c)
{
return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c);
return double3(a.x*c, a.y*c, a.z*c);
}
inline double operator*(const double3& a, const double3& b)
{
return a.data[0]*b.data[0]+a.data[1]*b.data[1]+a.data[2]*b.data[2];
return a.x*b.x+a.y*b.y+a.z*b.z;
}
inline double3 operator/(const double3& a, const double& c)
{
return double3(a.data[0]/c, a.data[1]/c, a.data[2]/c);
return double3(a.x/c, a.y/c, a.z/c);
}
inline double3 operator+(const double3& a, const double3& b)
{
return double3(a.data[0]+b.data[0], a.data[1]+b.data[1], a.data[2]+b.data[2]);
return double3(a.x+b.x, a.y+b.y, a.z+b.z);
}
inline double3 operator-(const double3& a, const double3& b)
{
return double3(a.data[0]-b.data[0], a.data[1]-b.data[1], a.data[2]-b.data[2]);
return double3(a.x-b.x, a.y-b.y, a.z-b.z);
}

6
io.cpp
View file

@ -174,9 +174,9 @@ void h5_write(const std::string file_name, const int step_num, const int N, cons
write_dataset("X", 2, (double*)x);
write_dataset("V", 2, (double*)v);
bool write_pot = (extra_mode ) % 2;
bool write_acc = (extra_mode >> 1) % 2;
bool write_jrk = (extra_mode >> 2) % 2;
bool write_pot = (extra_mode ) & 1;
bool write_acc = (extra_mode >> 1) & 1;
bool write_jrk = (extra_mode >> 2) & 1;
if (write_pot) write_dataset("POT", 1, (double*)pot);
if (write_acc) write_dataset("ACC", 2, (double*)acc);
if (write_jrk) write_dataset("JRK", 2, (double*)jrk);

424
phigrape.cpp
View file

@ -84,7 +84,6 @@ Config *config;
#include <numeric>
#include <stdexcept>
#define G6_NPIPE 2048
#include "grape6.h"
#include <mpi.h>
@ -165,6 +164,137 @@ void get_CPU_time(double *time_real, double *time_user, double *time_syst)
*time_user = *time_real;
}
void two_body_gravity(
const double m1, const double3 x1, const double3 v1,
const double m2, const double3 x2, const double3 v2,
const double eps,
double& pot1, double3& acc1, double3& jrk1,
double& pot2, double3& acc2, double3& jrk2)
{
double3 dx = x1 - x2;
double3 dv = v1 - v2;
double r2 = dx.norm2() + eps*eps;
double r = sqrt(r2);
double r3 = r2*r;
double r4 = r2*r2;
double RP = 3*(dx*dv)/r;
pot1 = -m2/r;
pot2 = -m1/r;
acc1 = -m2*dx/r3;
acc2 = m1*dx/r3;
jrk1 = -m2*(dv/r3 - RP*dx/r4);
jrk2 = m1*(dv/r3 - RP*dx/r4);
}
class Black_hole_physics {
public:
Black_hole_physics()
: count(0) {}
Black_hole_physics(const double m1, const double m2, const int myRank, const int rootRank)
: m1(m1), m2(m2), count(2), myRank(myRank), rootRank(rootRank) {}
void set_post_newtonian(const double c, const int pn_usage[7])
{
this->c = c;
std::copy(pn_usage, pn_usage+7, this->pn_usage);
}
void set_spins(const double spin1[3], const double spin2[3])
{
std::copy(spin1, spin1+3, this->spin1);
std::copy(spin2, spin2+3, this->spin2);
}
void adjust_softening(
const double eps_old, const double eps_new,
const double3& x1, const double3& x2,
const double3& v1, const double3& v2,
double& pot1, double& pot2,
double3& acc1, double3& acc2,
double3& jrk1, double3& jrk2);
void adjust_post_newtonian(
const double dt_bh, // pn_usage should be const
double3& x1, double3& x2,
double3& v1, double3& v2,
double& pot1, double& pot2,
double3& acc1, double3& acc2,
double3& jrk1, double3& jrk2);
public: //TODO make private
double m1, m2;
int count;
int myRank, rootRank;
double c;
int pn_usage[7];
double3 a_pn1[7], a_pn2[7], adot_pn1[7], adot_pn2[7];
double spin1[3], spin2[3];
};
void Black_hole_physics::adjust_softening(
const double eps_old, const double eps_new,
const double3& x1, const double3& x2,
const double3& v1, const double3& v2,
double& pot1, double& pot2,
double3& acc1, double3& acc2,
double3& jrk1, double3& jrk2)
{
double pot_bh1, pot_bh2;
double3 a_bh1, a_bh2, adot_bh1, adot_bh2;
// calculate and "minus" the BH <-> BH softened pot, acc & jerk
two_body_gravity(
m1, x1, v1,
m2, x2, v2,
eps_old,
pot_bh1, a_bh1, adot_bh1,
pot_bh2, a_bh2, adot_bh2);
pot1 -= pot_bh1;
pot2 -= pot_bh2;
acc1 -= a_bh1;
acc2 -= a_bh2;
jrk1 -= adot_bh1;
jrk2 -= adot_bh2;
// calculate and "plus" the new BH <-> BH unsoftened pot, acc, jerk
two_body_gravity(
m1, x1, v1,
m2, x2, v2,
eps_new,
pot_bh1, a_bh1, adot_bh1,
pot_bh2, a_bh2, adot_bh2);
pot1 += pot_bh1;
pot2 += pot_bh2;
acc1 += a_bh1;
acc2 += a_bh2;
jrk1 += adot_bh1;
jrk2 += adot_bh2;
}
void Black_hole_physics::adjust_post_newtonian(
const double dt_bh, // pn_usage should be const
double3& x1, double3& x2,
double3& v1, double3& v2,
double& pot1, double& pot2,
double3& acc1, double3& acc2,
double3& jrk1, double3& jrk2)
{
// calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
int tmp;
tmp = calc_force_pn_BH(m1, x1, v1, spin1,
m2, x2, v2, spin2,
c, dt_bh, pn_usage,
(double(*)[3])a_pn1, (double(*)[3])adot_pn1,
(double(*)[3])a_pn2, (double(*)[3])adot_pn2, myRank, rootRank);
if (tmp == 505) {
exit(-1); // Very ugly way to terminate
}
acc1 += a_pn1[1] + a_pn1[2] + a_pn1[3] + a_pn1[4] + a_pn1[5] + a_pn1[6];
acc2 += a_pn2[1] + a_pn2[2] + a_pn2[3] + a_pn2[4] + a_pn2[5] + a_pn2[6];
jrk1 += adot_pn1[1] + adot_pn1[2] + adot_pn1[3] + adot_pn1[4] + adot_pn1[5] + adot_pn1[6];
jrk2 += adot_pn2[1] + adot_pn2[2] + adot_pn2[3] + adot_pn2[4] + adot_pn2[5] + adot_pn2[6];
}
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) {
@ -679,59 +809,25 @@ inline void adjust_bsmbh_softening(const double eps, const double eps_bh, const
int main(int argc, char *argv[])
{
int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank,
diskstep=0, power, jjj, iii,
skip_con=0, tmp_i;
int skip_con=0; // skip_con should just be a static in the energy control function
double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
dt_bh, t_bh=0.0,
t_end, time_cur, dt_min, dt_max, min_t, min_t_loc,
eta_s, eta, eta_bh,
E_tot_0, E_tot_corr_0, E_tot_corr_sd_0,
double E_tot_0, E_tot_corr_0, E_tot_corr_sd_0,
rcm_sum=0.0, vcm_sum=0.0,
eps,
a2_mod, adot2_mod,
dt_tmp, dt2half, dt3over6,
timesteps=0.0, n_act_sum=0.0, n_act_distr[N_MAX];
double3 xcm, vcm, mom,
xdc, vdc,
a2, a3, a2dot1;
char processor_name[MPI_MAX_PROCESSOR_NAME],
inp_fname[30];
/* global variables */
Source_particle_list source_particle_list(N_MAX);
int N;
int *ind = source_particle_list.ind;
double *m = source_particle_list.m;
double3 *x = source_particle_list.x;
double3 *v = source_particle_list.v;
double *pot = source_particle_list.pot;
double3 *a = source_particle_list.a;
double3 *adot = source_particle_list.adot;
double *t = source_particle_list.t;
double *dt = source_particle_list.dt;
/* local variables */
int n_loc;
double3 xcm, vcm, xdc, vdc; // these should go away
/* data for active particles */
int n_act, ind_act[N_MAX];
double t_act[N_MAX],
pot_act_new[N_MAX],
pot_act_tmp[N_MAX];
double3 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;
// x_act_new and v_act_new arrays hold the predicted position and velocity of i-particles, which is later corrected before moving into the j-particle memory. The _old arrays hold the previous values which are needed to calculate the high derivatives. The _old arrays hold the previous values which are needed to calculate the high derivatives. The [a,adot]_act_tmp arrays hold the calculation results from each node. The [a,adot]_act_new arrays hold the reduced calculation results from all nodes.
double3 x_act_new[N_MAX], v_act_new[N_MAX],
a_act_old[N_MAX], adot_act_old[N_MAX],
a_act_tmp[N_MAX], adot_act_tmp[N_MAX],
a_act_new[N_MAX], adot_act_new[N_MAX];
double pot_ext[N_MAX], pot_act_ext[N_MAX]; // for all EXTPOT
@ -742,10 +838,6 @@ int main(int argc, char *argv[])
int inf_event[N_MAX];
double3 x_bbhc, v_bbhc;
/* some local settings for G6a board's */
int clusterid, numGPU, npipe=G6_NPIPE;
int new_tunit=51, new_xunit=51;
double3 zeros = {0, 0, 0}; // Dummy; can't really be const because of the GRAPE interface.
@ -756,13 +848,15 @@ int main(int argc, char *argv[])
/* Init MPI */
MPI_Init(&argc, &argv);
/* Define the total number of processors */
/* Define the total number of processors and the Rank of each processors */
int n_proc, myRank;
const int rootRank = 0;
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 */
int name_proc;
char processor_name[MPI_MAX_PROCESSOR_NAME];
MPI_Get_processor_name(processor_name, &name_proc);
/* Print the Rank and the names of processors */
@ -771,10 +865,23 @@ int main(int argc, char *argv[])
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* read the input parameters to the rootRank */
/* read the input parameters */
config = new Config("phigrape.conf");
Source_particle_list source_particle_list(N_MAX);
int *ind = source_particle_list.ind;
double *m = source_particle_list.m;
double3 *x = source_particle_list.x;
double3 *v = source_particle_list.v;
double *pot = source_particle_list.pot;
double3 *a = source_particle_list.a;
double3 *adot = source_particle_list.adot;
double *t = source_particle_list.t;
double *dt = source_particle_list.dt;
int diskstep, N;
double time_cur;
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");
@ -784,15 +891,14 @@ int main(int argc, char *argv[])
}
else
ascii_read(config->input_file_name, &diskstep, &N, &time_cur, m, x, v);
std::iota(ind, ind+N, 0);
eps = config->eps;
eta = config->eta;
t_end = config->t_end;
dt_disk = config->dt_disk;
dt_contr = config->dt_contr;
dt_bh = config->dt_bh;
double eps = config->eps;
double eta = config->eta;
double t_end = config->t_end;
double dt_disk = config->dt_disk;
double dt_contr = config->dt_contr;
double dt_bh = config->dt_bh;
if (myRank == rootRank) {
if (config->binary_smbh_influence_sphere_output) for (int i=0; i<N; i++) inf_event[i] = 0; // WARNING N wasn't set yet!
@ -807,18 +913,13 @@ int main(int argc, char *argv[])
printf("eta = %.6E \n", eta);
printf("\n");
printf("External Potential: \n");
printf("\n");
if ((diskstep == 0) && (time_cur == 0.0)) {
out = fopen("contr.dat", "w");
if ((diskstep == 0) && (time_cur == 0)) {
FILE *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);
@ -831,10 +932,7 @@ int main(int argc, char *argv[])
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);
@ -861,7 +959,15 @@ int main(int argc, char *argv[])
external_gravity_components.push_back(&ext_halo_plummer);
external_gravity_components.push_back(&ext_log_halo);
external_gravity_components.push_back(&ext_dehnen);
bool has_external_gravity = false;
for (auto component : external_gravity_components) {
if (component->is_active) {
has_external_gravity = true;
break;
}
}
if (has_external_gravity) printf("External Potential: \n\n");
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);
@ -870,22 +976,14 @@ int main(int argc, char *argv[])
printf("\n");
fflush(stdout);
eta_s = eta/ETA_S_CORR;
eta_bh = eta/ETA_BH_CORR;
double eta_bh = eta/ETA_BH_CORR;
dt_min = 1.0*pow(2.0, DTMINPOWER);
dt_max = 1.0*pow(2.0, DTMAXPOWER);
const double dt_min = pow(2, DTMINPOWER);
const double dt_max = std::min({dt_disk, dt_contr, pow(2, 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));
double t_disk = dt_disk*(1+floor(time_cur/dt_disk));
double t_contr = dt_contr*(1+floor(time_cur/dt_contr));
double t_bh = dt_bh*(1+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);
@ -896,15 +994,12 @@ int main(int argc, char *argv[])
std::fill(t, t+N, time_cur);
std::fill(dt, dt+N, dt_min);
n_loc = N/n_proc;
Calc_self_grav calc_self_grav(N, n_loc, clusterid, npipe, eps);
Active_search active_search(myRank, n_proc, n_loc, N);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
/* some local settings for G6a boards */
int clusterid, numGPU;
if (config->devices_per_node==0) {
MPI_Comm shmcomm;
MPI_Comm_split_type(MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, &shmcomm);
@ -919,10 +1014,23 @@ int main(int argc, char *argv[])
/* init the local GRAPEs */
g6_open(clusterid);
npipe = g6_npipes();
int npipe = g6_npipes();
g6_set_tunit(new_tunit);
g6_set_xunit(new_xunit);
int n_loc = N/n_proc;
Calc_self_grav calc_self_grav(N, n_loc, clusterid, npipe, eps);
Active_search active_search(myRank, n_proc, n_loc, N);
Black_hole_physics black_hole_physics;
if (config->live_smbh_count == 1)
black_hole_physics = Black_hole_physics(m[0], 0, myRank, rootRank);
else if (config->live_smbh_count == 2)
black_hole_physics = Black_hole_physics(m[0], m[1], myRank, rootRank);
if (config->binary_smbh_pn) {
black_hole_physics.set_post_newtonian(config->pn_c, config->pn_usage.data());
if (config->pn_usage[6]) black_hole_physics.set_spins(config->smbh1_spin.data(), config->smbh2_spin.data());
}
#ifdef ETICS
grapite_read_particle_tags(N, config->grapite_mask_file_name.c_str(), myRank, n_loc);
grapite_set_dt_exp(config->dt_scf);
@ -945,33 +1053,23 @@ int main(int argc, char *argv[])
grapite_cep_index = grapite_get_cep_index();
if (grapite_cep_index >= 0) {
grapite_calc_center(N, m, (double(*)[3])x, (double(*)[3])v, xcm, vcm, xdc, vdc);
memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
memcpy(v[grapite_cep_index], vdc, 3*sizeof(double));
x[grapite_cep_index] = xdc;
v[grapite_cep_index] = vdc;
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]; // isn't it just i?
iii = ind_act[i];
t_act[i] = t[iii];
} /* i */
calc_self_grav(time_cur, n_act, ind_act, x, v,
pot_act_tmp, a_act_tmp, adot_act_tmp);
calc_self_grav(time_cur, N, ind, x, v, pot_act_tmp, a_act_tmp, adot_act_tmp);
/* 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);
MPI_Allreduce(pot_act_tmp, pot, N, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(a_act_tmp, a, 3*N, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(adot_act_tmp, adot, 3*N, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
@ -980,57 +1078,13 @@ int main(int argc, char *argv[])
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
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
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->live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(eps, config->live_smbh_custom_eps, x[0], x[1], v[0], v[1], pot[0], pot[1], a[0], a[1], adot[0], adot[1]);
if (config->binary_smbh_pn) {
// calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
tmp_i = calc_force_pn_BH(m[i_bh1], x[i_bh1], v[i_bh1], config->smbh1_spin.data(),
int tmp;
tmp = calc_force_pn_BH(m[i_bh1], x[i_bh1], v[i_bh1], config->smbh1_spin.data(),
m[i_bh2], x[i_bh2], v[i_bh2], config->smbh2_spin.data(),
config->pn_c, dt[i_bh1], config->pn_usage.data(),
(double(*)[3])a_pn1, (double(*)[3])adot_pn1,
@ -1043,7 +1097,7 @@ int main(int argc, char *argv[])
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) {
if (tmp == 505) {
printf("PN RSDIST: %.8E \t %.8E \n", timesteps, time_cur);
fflush(stdout);
exit(-1);
@ -1052,7 +1106,7 @@ int main(int argc, char *argv[])
}
}
calc_ext_grav(external_gravity_components, n_act, x, v, pot_ext, a, adot);
calc_ext_grav(external_gravity_components, N, x, v, pot_ext, a, adot);
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
@ -1074,8 +1128,8 @@ int main(int argc, char *argv[])
if (grapite_cep_index >= 0) {
grapite_calc_center(N, m, (double(*)[3])x, (double(*)[3])v, xcm, vcm, xdc, vdc);
memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
memcpy(v[grapite_cep_index], vdc, 3*sizeof(double));
x[grapite_cep_index] = xdc;
v[grapite_cep_index] = vdc;
grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
}
#endif
@ -1086,13 +1140,14 @@ int main(int argc, char *argv[])
/* Define initial timestep for all particles on all nodes */
for (int j=0; j<N; j++) {
a2_mod = a[j].norm2();
adot2_mod = adot[j].norm2();
double a2_mod = a[j].norm2();
double adot2_mod = adot[j].norm2();
double dt_tmp, eta_s = eta/ETA_S_CORR;
if (adot2_mod==0) dt_tmp = eta_s; // That's weird, when will we have such a case?
else dt_tmp = eta_s*sqrt(a2_mod/adot2_mod);
power = log(dt_tmp)/log(2.0) - 1;
int power = log(dt_tmp)/log(2.0) - 1;
dt_tmp = pow(2.0, (double)power);
@ -1180,7 +1235,7 @@ int main(int argc, char *argv[])
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
#endif
min_t = active_search.get_minimum_time(t, dt);
double min_t = active_search.get_minimum_time(t, dt);
#ifdef TIMING
get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
@ -1247,8 +1302,8 @@ int main(int argc, char *argv[])
// NOTICE Why do we need pot_act and pot_act_ext? Probably redundant.
//pot_act[i] = pot[j_act];
pot_act_ext[i] = pot_ext[j_act]; // TODO almost certainly redundant
a_act[i] = a[j_act]; // TODO is used by calculate_high_derivatives, but should try to avoid!
adot_act[i] = adot[j_act]; // TODO is used by calculate_high_derivatives, but should try to avoid!
a_act_old[i] = a[j_act];
adot_act_old[i] = adot[j_act];
} /* i */
#ifdef TIMING
@ -1265,8 +1320,8 @@ int main(int argc, char *argv[])
for (int i=0; i<n_act; i++) {
int j_act = ind_act[i];
double dt = min_t - t[j_act];
dt2half = 0.5*dt*dt;
dt3over6 = (1./3.)*dt*dt2half;
double dt2half = 0.5*dt*dt;
double dt3over6 = (1./3.)*dt*dt2half;
x_act_new[i] = x[j_act] + v[j_act]*dt + a[j_act]*dt2half + adot[j_act]*dt3over6;
v_act_new[i] = v[j_act] + a[j_act]*dt + adot[j_act]*dt2half;
} /* i */
@ -1306,14 +1361,20 @@ int main(int argc, char *argv[])
#endif
if (config->live_smbh_count == 2) {
if (config->live_smbh_custom_eps >= 0) {
adjust_bsmbh_softening(eps, config->live_smbh_custom_eps, i_bh1, i_bh2, m[0], m[1], x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new);
}
//if (config->live_smbh_custom_eps >= 0) adjust_bsmbh_softening(eps, config->live_smbh_custom_eps, i_bh1, i_bh2, m[0], m[1], x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new);
if (config->live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(eps, config->live_smbh_custom_eps, x_act_new[i_bh1], x_act_new[i_bh2], v_act_new[i_bh1], v_act_new[i_bh2], pot_act_new[i_bh1], pot_act_new[i_bh2], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]);
//TODO the below works but it interferes with the printing
/* if (config->binary_smbh_pn) black_hole_physics.adjust_post_newtonian(
dt[i_bh1],
x_act_new[i_bh1], x_act_new[i_bh2],
v_act_new[i_bh1], v_act_new[i_bh2],
pot_act_new[i_bh1], pot_act_new[i_bh2], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]); */
if (config->binary_smbh_pn) {
// calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk
tmp_i = calc_force_pn_BH(m[0], x_act_new[i_bh1], v_act_new[i_bh1], config->smbh1_spin.data(),
int tmp;
tmp = calc_force_pn_BH(m[0], x_act_new[i_bh1], v_act_new[i_bh1], config->smbh1_spin.data(),
m[1], x_act_new[i_bh2], v_act_new[i_bh2], config->smbh2_spin.data(),
config->pn_c, dt[i_bh1], config->pn_usage.data(),
(double(*)[3])a_pn1, (double(*)[3])adot_pn1,
@ -1326,7 +1387,7 @@ int main(int argc, char *argv[])
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) {
if (tmp == 505) {
printf("PN RSDIST: TS = %.8E \t t = %.8E \n", timesteps, time_cur);
fflush(stdout);
exit(-1);
@ -1348,10 +1409,10 @@ int main(int argc, char *argv[])
// NOTICE looks like we're doing three unrelated things in this loop: (1) correcting positions and velocities (2) calculating new steps, and (3) putting the corrected values from the _act_new back in the _act arrays.
// After going back to the _act arrays they don't do much before they go back to the main arrays, so this copy seems redundant (the SMBH influence sphere printout needs these values but it should be a function anyway).
// TODO split this loop into the three tasks it is doing, and remove the redundancy.
dt_tmp = min_t - t_act[i];
double dt_tmp = min_t - t_act[i];
double3 a2, a3;
calc_high_derivatives(dt_tmp, a_act[i], a_act_new[i], adot_act[i], adot_act_new[i], &a2, &a3);
calc_high_derivatives(dt_tmp, a_act_old[i], a_act_new[i], adot_act_old[i], adot_act_new[i], &a2, &a3);
corrector(dt_tmp, a2, a3, &x_act_new[i], &v_act_new[i]);
@ -1365,7 +1426,7 @@ int main(int argc, char *argv[])
//TODO the below should be moved to a function
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;
int power = log(dt_new)/log(2.0) - 1;
dt_tmp = pow(2.0, (double)power); // TODO why is this casting needed here?
}
@ -1420,13 +1481,12 @@ int main(int argc, char *argv[])
#ifdef ETICS
if (ind_act[i] == grapite_cep_index) grapite_update_cep(t[grapite_cep_index], x[grapite_cep_index], v[grapite_cep_index], a[grapite_cep_index], adot[grapite_cep_index]); // All ranks should do it.
#endif
cur_rank = ind_act[i]/n_loc;
int cur_rank = ind_act[i]/n_loc;
if (myRank == cur_rank) {
int j_act = ind_act[i];
jjj = ind_act[i] - myRank*n_loc;
g6_set_j_particle(clusterid, jjj, ind_act[i], t[j_act], dt[j_act], m[j_act], zeros, adot[j_act]*(1./6.), a[j_act]*0.5, v[j_act], x[j_act]);
int address = ind_act[i] - myRank*n_loc;
g6_set_j_particle(clusterid, address, ind_act[i], t[j_act], dt[j_act], m[j_act], zeros, adot[j_act]*(1./6.), a[j_act]*0.5, v[j_act], x[j_act]);
} /* if (myRank == cur_rank) */
@ -1471,7 +1531,7 @@ int main(int argc, char *argv[])
/* possible OUT for timing !!! */
#ifdef TIMING
out = fopen("timing.dat", "a");
FILE *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 +
@ -1500,8 +1560,8 @@ int main(int argc, char *argv[])
// 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.
if (grapite_cep_index >= 0) {
grapite_calc_center(N, m, (double(*)[3])x, (double(*)[3])v, xcm, vcm, xdc, vdc);
memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
memcpy(v[grapite_cep_index], vdc, 3*sizeof(double));
x[grapite_cep_index] = xdc;
v[grapite_cep_index] = vdc;
grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
}
#endif