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

Started moving toward dynamical allocation; removed the _loc variables which were only used in the zeroth step and were pretty unnecessary anyway

Browse source
This commit is contained in:
Yohai Meiron 2020-04-03 21:40:11 -04:00
parent ea94dbb626
commit 99ec3064ed
2 changed files with 59 additions and 118 deletions
Download patch file Download diff file Expand all files Collapse all files

2
config.cpp
View file

@ -5,6 +5,8 @@
#include <algorithm> #include <algorithm>
#include <limits> #include <limits>
// Would be a bit more elegant to do the whole thing with std::variant.
using Dictionary = std::unordered_map<std::string,std::string>; using Dictionary = std::unordered_map<std::string,std::string>;
static constexpr double nix = -std::numeric_limits<double>::max(); // avoid nans static constexpr double nix = -std::numeric_limits<double>::max(); // avoid nans

173
phigrape.cpp
View file

@ -77,25 +77,6 @@ Config *config;
#define DTMAXPOWER -3.0 #define DTMAXPOWER -3.0
#define DTMINPOWER -36.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 <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <math.h> #include <math.h>
@ -114,46 +95,19 @@ Config *config;
#include <mpi.h> #include <mpi.h>
/* Some "good" functions and constants... */ /* Some "good" functions and constants... */
#define SIG(x) (((x)<0) ? (-1):(1) ) // TODO replace with inline functions
#define ABS(x) (((x)<0) ? (-x):(x) )
#define MAX(a,b) (((a)>(b)) ? (a):(b) )
#define MIN(a,b) (((a)<(b)) ? (a):(b) ) #define MIN(a,b) (((a)<(b)) ? (a):(b) )
#define SQR(x) ((x)*(x) ) #define SQR(x) ((x)*(x) )
#define POW3(x) ((x)*SQR(x) )
#define Pi 3.14159265358979323846
#define TWOPi 6.283185307179
#define sqrt_TWOPi 2.506628274631
/*
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 KB 1024
#define MB (KB*KB) #define MB (KB*KB)
#define N_MAX (6*MB) #define N_MAX (6*MB)
#define N_MAX_loc (2*MB)
double L[3]; // needed in pn_bh_spin.c double L[3]; // needed in pn_bh_spin.c
// Needed for things related to BHs // Needed for things related to BHs
#include "debug.h" #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_real0, CPU_time_user0, CPU_time_syst0;
double CPU_time_real, CPU_time_user, CPU_time_syst; double CPU_time_real, CPU_time_user, CPU_time_syst;
@ -185,11 +139,6 @@ double3 x_i[G6_NPIPE], v_i[G6_NPIPE],
int new_tunit=51, new_xunit=51; int new_tunit=51, new_xunit=51;
double ti=0.0; double ti=0.0;
double3 a2by18, a1by6, aby2;
/* normalization... */
double m_norm, r_norm, v_norm, t_norm;
double eps_BH=0.0; double eps_BH=0.0;
@ -200,14 +149,10 @@ double3 x_bh1, x_bh2, v_bh1, v_bh2;
double pot_bh1, pot_bh2; double pot_bh1, pot_bh2;
double3 a_bh1, adot_bh1, a_bh2, adot_bh2; double3 a_bh1, adot_bh1, a_bh2, adot_bh2;
//double eps_BH = 0.0;
#include "n_bh.c" #include "n_bh.c"
double3 a_pn1[7], adot_pn1[7], a_pn2[7], adot_pn2[7]; double3 a_pn1[7], adot_pn1[7], a_pn2[7], adot_pn2[7];
#include "pn_bh_spin.c" #include "pn_bh_spin.c"
#ifdef ETICS #ifdef ETICS
@ -318,6 +263,7 @@ void write_bh_nb_data(double time_cur, int N, double m[], double3 x[], double3 v
int nb = config->live_smbh_neighbor_number; int nb = config->live_smbh_neighbor_number;
int ind_sort[N_MAX]; int ind_sort[N_MAX];
double var_sort[N_MAX]; double var_sort[N_MAX];
//TODO you don't want and probably don't need to allocate these here. Maybe just need size nb, or maybe use as private variables.
auto out = fopen("bh_neighbors.dat", "a"); auto out = fopen("bh_neighbors.dat", "a");
@ -583,6 +529,39 @@ private:
int *ind_act_loc; int *ind_act_loc;
}; };
class Source_particle_list {
public:
Source_particle_list(int N)
: N(N)
{
ind = new int[N];
m = new double[N];
x = new double3[N];
v = new double3[N];
pot = new double[N];
a = new double3[N];
adot = new double3[N];
t = new double[N];
dt = new double[N];
}
~Source_particle_list()
{
delete[] ind;
delete[] m;
delete[] x;
delete[] v;
delete[] pot;
delete[] a;
delete[] adot;
delete[] t;
delete[] dt;
}
int N;
int *ind;
double3 *x, *v, *a, *adot;
double *m, *t, *dt, *pot;
};
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank, int name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank,
@ -612,16 +591,21 @@ int main(int argc, char *argv[])
/* global variables */ /* global variables */
int N, ind[N_MAX]; Source_particle_list source_particle_list(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]; 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 */ /* local variables */
int n_loc;
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 */ /* data for active particles */
@ -656,6 +640,7 @@ int main(int argc, char *argv[])
double s_bh1[3] = {0.0, 0.0, 1.0}; double s_bh1[3] = {0.0, 0.0, 1.0};
double s_bh2[3] = {0.0, 0.0, 1.0}; double s_bh2[3] = {0.0, 0.0, 1.0};
double3 zeros = {0, 0, 0}; // Dummy; can't really be const because of the GRAPE interface.
/* INIT the rand() !!! */ /* INIT the rand() !!! */
srand(19640916); /* it is just my birthday :-) */ srand(19640916); /* it is just my birthday :-) */
@ -848,14 +833,6 @@ int main(int argc, char *argv[])
/* Wait to all processors to finish his works... */ /* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD); 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... */ /* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
@ -885,17 +862,11 @@ int main(int argc, char *argv[])
grapite_set_t_exp(t_exp); grapite_set_t_exp(t_exp);
#endif #endif
/* initial load the particles to the local GRAPE's */
nj=n_loc;
/* load the nj particles to the G6 */ /* load the nj particles to the G6 */
for (int j=0; j<nj; j++) { for (int j=0; j<n_loc; j++) {
a2by18 = {0, 0, 0}; int kkk = j + myRank*n_loc;
a1by6 = {0, 0, 0}; g6_set_j_particle(clusterid, j, ind[kkk], t[kkk], dt[kkk], m[kkk], zeros, zeros, zeros, v[kkk], x[kkk]);
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 */ } /* j */
#ifdef ETICS #ifdef ETICS
@ -916,12 +887,6 @@ int main(int argc, char *argv[])
memcpy(x[grapite_cep_index], xdc, 3*sizeof(double)); memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
memcpy(v[grapite_cep_index], vdc, 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! double zeros[3] = {0., 0., 0.}; // We haven't calculated the force yet!
grapite_update_cep(time_cur, xdc, vdc, zeros, zeros); grapite_update_cep(time_cur, xdc, vdc, zeros, zeros);
#endif #endif
@ -1064,11 +1029,6 @@ int main(int argc, char *argv[])
/* Wait to all processors to finish his works... */ /* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD); 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 #ifdef ETICS_CEP
// First calculate the DC // First calculate the DC
grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc); grapite_calc_center(N, m, x, v, xcm, vcm, xdc, vdc);
@ -1077,12 +1037,6 @@ int main(int argc, char *argv[])
memcpy(x[grapite_cep_index], xdc, 3*sizeof(double)); memcpy(x[grapite_cep_index], xdc, 3*sizeof(double));
memcpy(v[grapite_cep_index], vdc, 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]); grapite_update_cep(time_cur, xdc, vdc, a[grapite_cep_index], adot[grapite_cep_index]);
#endif #endif
@ -1125,23 +1079,12 @@ int main(int argc, char *argv[])
/* Wait to all processors to finish his works... */ /* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD); 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 */ /* load the new values for particles to the local GRAPE's */
nj=n_loc;
/* load the nj particles to the G6 */ /* load the nj particles to the G6 */
for (int j=0; j<nj; j++) { for (int j=0; j<n_loc; j++) {
a2by18 = {0, 0, 0}; int kkk = j + myRank*n_loc;
a1by6 = adot_loc[j]*(1./6.); g6_set_j_particle(clusterid, j, ind[kkk], t[kkk], dt[kkk], m[kkk], zeros, adot[kkk]*(1./6.), a[kkk]*0.5, v[kkk], x[kkk]);
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 */ } /* j */
if (myRank == rootRank) { if (myRank == rootRank) {
@ -1592,11 +1535,7 @@ int main(int argc, char *argv[])
jjj = ind_act[i] - myRank*n_loc; jjj = ind_act[i] - myRank*n_loc;
a2by18 = {0, 0, 0}; g6_set_j_particle(clusterid, jjj, ind_act[i], t_act[i], dt_act[i], m_act[i], zeros, adot_act[i]*(1./6.), a_act[i]*0.5, v_act[i], x_act[i]);
a1by6 = adot_act[i]*(1./6.);
aby2 = a_act[i]*0.5;
g6_set_j_particle(clusterid, jjj, ind_act[i], t_act[i], dt_act[i], m_act[i], a2by18, a1by6, aby2, v_act[i], x_act[i]);
} /* if (myRank == cur_rank) */ } /* if (myRank == cur_rank) */