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

Moved config to inside main and changed several arrays into vectors

Browse source
This commit is contained in:
Yohai Meiron 2020-10-27 23:55:55 -04:00
parent 329dd2ca4d
commit 747f9f9d89
7 changed files with 115 additions and 108 deletions
Download patch file Download diff file Expand all files Collapse all files

167
phigrape.cpp
View file

@ -28,7 +28,6 @@
#include "grapite.h"
#endif
Config *config;
//chrono::steady_clock::time_point walltime_start;
namespace std::chrono {
@ -58,7 +57,7 @@ public:
acc_loc.resize(N);
jrk_loc.resize(N);
}
void operator()(const double t, const int n_act, std::vector<int>& ind_act, const double3 x_act[], const double3 v_act[],
void operator()(const double t, const int n_act, std::vector<int> &ind_act, std::vector<double3> &x_act, std::vector<double3> &v_act,
std::vector<double>& pot, double3 acc[], double3 jrk[])
{
g6_set_ti(clusterid, t);
@ -84,7 +83,7 @@ private:
std::vector<double3> acc_loc, jrk_loc;
};
void calc_ext_grav(std::vector<External_gravity*> &external_gravity_components, int n, double3 *x, double3 *v, double *pot, double3 *acc, double3* jrk)
void calc_ext_grav(std::vector<External_gravity*> &external_gravity_components, int n, const std::vector<double3> &x, const std::vector<double3> &v, double *pot, double3 *acc, double3* jrk)
// TODO should just be a class that has this pointer array as a member
{
std::fill(pot, pot+n, 0.);
@ -94,7 +93,7 @@ void calc_ext_grav(std::vector<External_gravity*> &external_gravity_components,
}
}
void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, int N, double m[], double3 x[], double3 v[], const std::vector<double>& pot, double pot_ext[])
void energy_contr(const double time_cur, const double timesteps, const double n_act_sum, const double g6_calls, int N, const std::vector<double> &m, const std::vector<double3> &x, const std::vector<double3> &v, const std::vector<double> &pot, double pot_ext[])
{
double E_pot = 0;
for (int i=0; i<N; i++) E_pot += m[i]*pot[i];
@ -160,8 +159,8 @@ void energy_contr(const double time_cur, const double timesteps, const double n_
class Active_search {
// TODO you can add pointers to t and dt at the constructor, no point giving them at get_minimum_time but without the size.
public:
Active_search(const int myRank, const int n_proc, const int n_loc, const int N)
: myRank(myRank), n_proc(n_proc), n_loc(n_loc), N(N)
Active_search(const int myRank, const int n_proc, const int n_loc, const int N, bool grapite_active_search_flag)
: myRank(myRank), n_proc(n_proc), n_loc(n_loc), N(N), grapite_active_search_flag(grapite_active_search_flag)
{
ind_act_loc = new int[n_loc];
}
@ -170,7 +169,7 @@ public:
{
double min_t_loc, min_t;
#ifdef ETICS
if (config->grapite_active_search) {
if (grapite_active_search_flag) {
min_t_loc = grapite_get_minimum_time();
} else
#endif
@ -188,7 +187,7 @@ public:
void get_active_indices(const double min_t, const double t[], const double dt[], int ind_act[], int& n_act)
{
#ifdef ETICS
if (config->grapite_active_search) {
if (grapite_active_search_flag) {
int n_act_loc;
grapite_active_search(min_t, ind_act_loc, &n_act_loc);
if (myRank > 0)
@ -215,6 +214,7 @@ public:
private:
int myRank, n_proc, n_loc, N;
int *ind_act_loc;
bool grapite_active_search_flag;
};
inline void calc_high_derivatives(const double dt_tmp, const double3 a_old, const double3 a_new, const double3 a1_old, const double3 a1_new, double3& a2, double3& a3)
@ -254,6 +254,9 @@ inline double aarseth_step(const double eta, const double dt, const double3 a, c
int main(int argc, char *argv[])
{
/* read the input parameters */
const Config config("phigrape.conf");
timer.start();
double timesteps=0.0, n_act_sum=0.0;
@ -280,23 +283,21 @@ int main(int argc, char *argv[])
/* Print the Rank and the names of processors */
printf("Rank of the processor %03d on %s \n", myRank, processor_name);
/* read the input parameters */
config = new Config("phigrape.conf");
int diskstep, N;
double time_cur;
// The memory for m, x, and v is allocated inside h5_read or ascii_read
double *m;
double3 *x, *v;
if (is_hdf5(config->input_file_name)) {
//double *m;
std::vector<double> m;
std::vector<double3> x, v;
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);
h5_read(config.input_file_name, &diskstep, &N, &time_cur, m.data(), x.data(), v.data());
}
else
ascii_read(config->input_file_name, diskstep, N, time_cur, &m, &x, &v);
ascii_read(config.input_file_name, diskstep, N, time_cur, m, x, v);
std::vector<int> ind(N);
std::iota(begin(ind), end(ind), 0);
@ -308,16 +309,16 @@ int main(int argc, char *argv[])
int n_act;
std::vector<int> ind_act(N);
std::vector<double> pot_act_new(N);
std::vector<double3> x_act_new(N), v_act_new(N);
double *pot_act_ext = new double[N];
double3 *x_act_new = new double3[N], *v_act_new = new double3[N],
*a_act_new = new double3[N], *adot_act_new = new double3[N];
double3 *a_act_new = new double3[N], *adot_act_new = new double3[N];
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;
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) {
printf("\n");
@ -333,11 +334,11 @@ int main(int argc, char *argv[])
if ((diskstep == 0) && (time_cur == 0)) {
FILE *out = fopen("contr.dat", "w");
fclose(out);
if (config->live_smbh_output && (config->live_smbh_count > 0)) {
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)) {
if ((config.live_smbh_neighbor_output) && (config.live_smbh_count > 0)) {
out = fopen("bh_neighbors.dat", "w");
fclose(out);
}
@ -345,16 +346,16 @@ int main(int argc, char *argv[])
} /* if (myRank == rootRank) */
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);
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);
Dehnen ext_dehnen(config->ext_dehnen_m*normalization_mass, config->ext_dehnen_r*normalization_length, config->ext_dehnen_gamma);
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);
Dehnen ext_dehnen(config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma);
std::vector<External_gravity*> external_gravity_components;
external_gravity_components.push_back(&ext_bulge);
@ -371,11 +372,11 @@ int main(int argc, char *argv[])
}
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);
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);
if (ext_dehnen.is_active) printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n", config->ext_dehnen_m*normalization_mass, config->ext_dehnen_r*normalization_length, config->ext_dehnen_gamma);
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);
if (ext_dehnen.is_active) printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n", config.ext_dehnen_m*normalization_mass, config.ext_dehnen_r*normalization_length, config.ext_dehnen_gamma);
printf("\n");
fflush(stdout);
@ -399,13 +400,13 @@ int main(int argc, char *argv[])
/* some local settings for G6a boards */
int clusterid, numGPU;
if (config->devices_per_node==0) {
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);
MPI_Comm_size(shmcomm, &numGPU);
MPI_Comm_rank(shmcomm, &clusterid);
} else {
numGPU = config->devices_per_node;
numGPU = config.devices_per_node;
clusterid = myRank % numGPU;
}
printf("Rank of the processor %03d : Number of GPUs %01d : Cluster ID %01d \n", myRank, numGPU, clusterid);
@ -418,30 +419,33 @@ int main(int argc, char *argv[])
g6_set_xunit(51);
#ifdef ETICS
grapite_set_dev_exec_threshold(config->grapite_dev_exec_threshold);
grapite_set_dev_exec_threshold(config.grapite_dev_exec_threshold);
bool grapite_active_search_flag = config.grapite_active_search;
#else
bool grapite_active_search_flag = false;
#endif
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);
Active_search active_search(myRank, n_proc, n_loc, N, grapite_active_search_flag);
Black_hole_physics black_hole_physics;
if (config->live_smbh_count == 1)
if (config.live_smbh_count == 1)
black_hole_physics = Black_hole_physics(m[0], 0, myRank, rootRank);
else if (config->live_smbh_count == 2)
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());
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());
}
black_hole_physics.set_softening(eps, config->live_smbh_custom_eps);
black_hole_physics.set_softening(eps, config.live_smbh_custom_eps);
Binary_smbh_influence_sphere_output binary_smbh_influence_sphere_output(config->binary_smbh_influence_radius_factor, N, m, x, v, pot, dt);
Binary_smbh_influence_sphere_output binary_smbh_influence_sphere_output(config.binary_smbh_influence_radius_factor, N, m, x, v, pot, dt);
Write_bh_nb_data write_bh_nb_data(config->live_smbh_neighbor_number, config->live_smbh_count, N, m, x, v);
Write_bh_nb_data write_bh_nb_data(config.live_smbh_neighbor_number, config.live_smbh_count, N, m, x, v);
#ifdef ETICS
grapite_read_particle_tags(N, config->grapite_mask_file_name.c_str(), myRank, n_loc);
grapite_set_dt_exp(config->dt_scf);
grapite_read_particle_tags(N, config.grapite_mask_file_name.c_str(), myRank, n_loc);
grapite_set_dt_exp(config.dt_scf);
grapite_set_t_exp(time_cur);
#endif
@ -470,10 +474,10 @@ int main(int argc, char *argv[])
/* define the all particles as a active on all the processors for the first time grav calc. */
calc_self_grav(time_cur, N, ind, x, v, pot, a, adot);
if (config->live_smbh_count == 2) {
if (config.live_smbh_count == 2) {
black_hole_physics.set_xv(x[0], x[1], v[0], v[1]);
if (config->live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(pot[0], pot[1], a[0], a[1], adot[0], adot[1]);
if (config->binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[0], a[0], a[1], adot[0], adot[1]);
if (config.live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(pot[0], pot[1], a[0], a[1], adot[0], adot[1]);
if (config.binary_smbh_pn) black_hole_physics.adjust_post_newtonian(dt[0], a[0], a[1], adot[0], adot[1]);
}
calc_ext_grav(external_gravity_components, N, x, v, pot_ext, a, adot);
@ -482,7 +486,7 @@ int main(int argc, char *argv[])
if (myRank == rootRank) energy_contr(time_cur, timesteps, n_act_sum, calc_self_grav.g6_calls, N, m, x, v, pot, pot_ext);
#ifdef ETICS
if (config->etics_dump_coeffs && (diskstep==0)) {
if (config.etics_dump_coeffs && (diskstep==0)) {
char out_fname[256];
sprintf(out_fname, "coeffs.%06d.%02d.dat", 0, myRank);
grapite_dump(out_fname, 2);
@ -515,7 +519,7 @@ int main(int argc, char *argv[])
dt[j] = dt_tmp;
if (config->dt_min_warning && (myRank == 0)) {
if (config.dt_min_warning && (myRank == 0)) {
if (dt[j] == dt_min) {
printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[j], ind[j]);
fflush(stdout);
@ -523,9 +527,9 @@ int main(int argc, char *argv[])
}
} /* j */
if (config->live_smbh_count > 0) {
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;
for (int i=0; i<config.live_smbh_count; i++) dt[i] = min_dt;
}
/* load the new values for particles to the local GRAPEs */
@ -536,10 +540,10 @@ int main(int argc, char *argv[])
if (myRank == rootRank) {
/* Write BH data... */
if (config->live_smbh_output) black_hole_physics.write_bh_data(time_cur, m, x, v, pot, a, adot, dt);
if (config.live_smbh_output) black_hole_physics.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);
if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur);
} /* if (myRank == rootRank) */
@ -558,8 +562,8 @@ int main(int argc, char *argv[])
/* Find the BH(s) indices in the active list */
int i_bh1=0, i_bh2=1;
#ifdef ETICS
int n_bh = config->live_smbh_count;
if (config->grapite_active_search && (n_bh>0)) {
int n_bh = config.live_smbh_count;
if (config.grapite_active_search && (n_bh>0)) {
int act_def_grapite_bh_count = 0;
int i_bh[n_bh];
for (int i=0; i<n_act; i++) {
@ -586,10 +590,10 @@ int main(int argc, char *argv[])
/* Calculate gravity on active particles */
calc_self_grav(min_t, n_act, ind_act, x_act_new, v_act_new, pot_act_new, a_act_new, adot_act_new);
if (config->live_smbh_count == 2) {
if (config.live_smbh_count == 2) {
black_hole_physics.set_xv(x_act_new[i_bh1], x_act_new[i_bh2], v_act_new[i_bh1], v_act_new[i_bh2]);
if (config->live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(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) black_hole_physics.adjust_post_newtonian(dt[i_bh1], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]);
if (config.live_smbh_custom_eps >= 0) black_hole_physics.adjust_softening(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) black_hole_physics.adjust_post_newtonian(dt[i_bh1], a_act_new[i_bh1], a_act_new[i_bh2], adot_act_new[i_bh1], adot_act_new[i_bh2]);
}
/* Calculate gravity on active particles due to external forces */
@ -608,7 +612,7 @@ int main(int argc, char *argv[])
//TODO make beautiful
double eta_curr;
if ((config->live_smbh_count > 0) && (ind_act[i] < config->live_smbh_count)) eta_curr = eta_bh;
if ((config.live_smbh_count > 0) && (ind_act[i] < config.live_smbh_count)) eta_curr = eta_bh;
else eta_curr = eta;
double dt_new = aarseth_step(eta_curr, dt_tmp, a_act_new[i], adot_act_new[i], a2, a3);
@ -624,7 +628,7 @@ int main(int argc, char *argv[])
dt_tmp *= 2;
}
if (config->dt_min_warning && (myRank == 0)) {
if (config.dt_min_warning && (myRank == 0)) {
if (dt_tmp == dt_min) {
printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt_tmp, ind_act[i]);
fflush(stdout);
@ -643,12 +647,12 @@ int main(int argc, char *argv[])
} /* i */
/* define the min. dt over all the act. part. and set it also for the BH... */
if (config->live_smbh_count > 0) {
if (config->live_smbh_count>=1) dt[0] = min_dt;
if (config->live_smbh_count==2) dt[1] = min_dt;
if (config.live_smbh_count > 0) {
if (config.live_smbh_count>=1) dt[0] = min_dt;
if (config.live_smbh_count==2) dt[1] = min_dt;
}
if (config->binary_smbh_influence_sphere_output && (myRank == rootRank)) {
if (config.binary_smbh_influence_sphere_output && (myRank == rootRank)) {
//TODO clean up this mass. We don't want to have all these _act arrays; they are only needed for this lame function.
binary_smbh_influence_sphere_output(ind_act, n_act, timesteps, time_cur);
}
@ -675,10 +679,10 @@ int main(int argc, char *argv[])
if (time_cur >= t_bh) {
if (myRank == rootRank) {
/* Write BH data... */
if (config->live_smbh_output) black_hole_physics.write_bh_data(time_cur, m, x, v, pot, a, adot, dt);
if (config.live_smbh_output) black_hole_physics.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);
if (config.live_smbh_neighbor_output) write_bh_nb_data(time_cur);
} /* if (myRank == rootRank) */
@ -689,7 +693,7 @@ int main(int argc, char *argv[])
if (myRank == rootRank) {
energy_contr(time_cur, timesteps, n_act_sum, calc_self_grav.g6_calls, 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);
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);
} /* if (myRank == rootRank) */
@ -717,12 +721,12 @@ int main(int argc, char *argv[])
diskstep++;
if (myRank == rootRank) {
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, config->output_extra_mode, config->output_hdf5_double_precision);
else ascii_write(std::string(out_fname) + ".dat", diskstep, N, time_cur, m, x, v, config->output_ascii_precision);
if (config.output_hdf5) h5_write(std::string(out_fname) + ".h5", diskstep, N, time_cur, m, x, v, pot, a, adot, config.output_extra_mode, config.output_hdf5_double_precision);
else ascii_write(std::string(out_fname) + ".dat", diskstep, N, time_cur, m, x, v, config.output_ascii_precision);
} /* if (myRank == rootRank) */
#ifdef ETICS
if (config->etics_dump_coeffs) {
if (config.etics_dump_coeffs) {
sprintf(out_fname, "coeffs.%06d.%02d.dat", diskstep, myRank);
grapite_dump(out_fname, 2);
}
@ -746,8 +750,7 @@ int main(int argc, char *argv[])
fflush(stdout);
} /* if (myRank == rootRank) */
delete config;
delete[] m; delete[] x; delete[] v; delete[] a; delete[] adot; delete[] pot_ext; delete[] t; delete[] dt; delete[] x_act_new; delete[] v_act_new; delete[] a_act_new; delete[] adot_act_new; delete[] pot_act_ext;
delete[] a; delete[] adot; delete[] pot_ext; delete[] t; delete[] dt; delete[] a_act_new; delete[] adot_act_new; delete[] pot_act_ext;
/* Finalize the MPI work */
MPI_Finalize();