diff --git a/config.cpp b/config.cpp index 28ac76f..e9b001d 100644 --- a/config.cpp +++ b/config.cpp @@ -5,6 +5,8 @@ #include #include +// Would be a bit more elegant to do the whole thing with std::variant. + using Dictionary = std::unordered_map; static constexpr double nix = -std::numeric_limits::max(); // avoid nans diff --git a/phigrape.cpp b/phigrape.cpp index 0e61467..681089e 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -77,25 +77,6 @@ Config *config; #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 #include #include @@ -114,46 +95,19 @@ Config *config; #include /* 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) ) +// TODO replace with inline functions #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 - -/* - 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; @@ -185,11 +139,6 @@ double3 x_i[G6_NPIPE], v_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; @@ -200,14 +149,10 @@ double3 x_bh1, x_bh2, v_bh1, v_bh2; 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 @@ -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 ind_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"); @@ -583,6 +529,39 @@ private: 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 name_proc, n_proc=1, myRank=0, rootRank=0, cur_rank, @@ -612,16 +591,21 @@ int main(int argc, char *argv[]) /* 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]; + 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, 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]; + int n_loc; /* data for active particles */ @@ -655,7 +639,8 @@ int main(int argc, char *argv[]) double s_bh1[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() !!! */ 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... */ 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); @@ -885,17 +862,11 @@ int main(int argc, char *argv[]) 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