Moved more variables from global to main and removed some unused variables

phigrape.cpp

@@ -314,46 +314,6 @@ double L[3]; // needed in pn_bh_spin.c
 // double var_sort[N_MAX];
 
 
-double3 xcm, vcm, mom,
-        xdc, vdc,
-        a2, a3, a2dot1;
-
-char   processor_name[MPI_MAX_PROCESSOR_NAME],
-    inp_fname[30],
-    out_fname[30];
-
-/* global variables */
-
-int    N, N_star, N_bh,
-    ind[N_MAX], name[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],
-       pot_act_tmp_loc[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],
-        a_act_tmp_loc[N_MAX], adot_act_tmp_loc[N_MAX];
-
-FILE   *inp, *out;
-
 double CPU_time_real0, CPU_time_user0, CPU_time_syst0;
 double CPU_time_real,  CPU_time_user,  CPU_time_syst;
 
@@ -584,12 +544,11 @@ gettimeofday(&tv, NULL);
 
 }
 
-void read_data(int *diskstep, int *N, double *time_cur, int ind[], double m[], double3 x[], double3 v[])
+void read_data(char inp_fname[], int *diskstep, int *N, double *time_cur, int ind[], double m[], double3 x[], double3 v[])
 {
     auto inp = fopen(inp_fname, "r");
     fscanf(inp, "%d \n", diskstep);
     fscanf(inp, "%d \n", N);
-    printf("zzzzzz reading *%s* %d\n", inp_fname, *N);
     fscanf(inp, "%lE \n", time_cur);
     for (int i=0; i<*N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2]);
     fclose(inp);
@@ -597,7 +556,7 @@ void read_data(int *diskstep, int *N, double *time_cur, int ind[], double m[], d
 
 void write_snap_data(char out_fname[], int diskstep, int N, double time_cur, int ind[], double m[], double3 x[], double3 v[])
 {
-    out = fopen(out_fname, "w");
+    auto out = fopen(out_fname, "w");
     fprintf(out,"%06d \n", diskstep);
     fprintf(out,"%07d \n", N);
     fprintf(out,"%.10E \n", time_cur);
@@ -611,28 +570,7 @@ void write_snap_data(char out_fname[], int diskstep, int N, double time_cur, int
     fclose(out);
 }
 
-// // // // // // // // // // // void write_cont_data() // virtually identical to write_snap_data()
+void write_bh_data(double time_cur, double m[], double3 x[], double3 v[], double pot[], double3 a[], double3 adot[], double dt[])
-// // // // // // // // // // // {
-// // // // // // // // // // // out = fopen("data.con","w");
-// // // // // // // // // // // fprintf(out,"%06d \n", diskstep);
-// // // // // // // // // // // 
-// // // // // // // // // // // fprintf(out,"%07d \n", N);
-// // // // // // // // // // // 
-// // // // // // // // // // // fprintf(out,"%.16E \n", time_cur);
-// // // // // // // // // // // 
-// // // // // // // // // // // for (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[])
 {
     if (config->live_smbh_count == 2) {
 
@@ -735,7 +673,7 @@ void write_bh_data(double time_cur, double m[], double3 x[], double3 v[])
         fprintf(out,"\n");
         fclose(out);
     } else if (config->live_smbh_count == 1) {
-        out = fopen("bh.dat", "a");
+        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]) );
@@ -753,13 +691,13 @@ void write_bh_data(double time_cur, double m[], double3 x[], double3 v[])
     }
 }
 
-void write_bh_nb_data(double time_cur, double m[], double3 x[], double3 v[])
+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];
 
-    out = fopen("bh_neighbors.dat", "a");
+    auto out = fopen("bh_neighbors.dat", "a");
 
     /* 1st BH */
 
@@ -1039,7 +977,13 @@ for (i=0; i<ni; i++)
 
 }
 
-void calc_self_grav(double t, double eps2, double &g6_calls)
+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 */
 
@@ -1047,6 +991,11 @@ void calc_self_grav(double t, double eps2, double &g6_calls)
     get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
 #endif
 
+// // // // // // //     /* bad variables... */
+// // // // // // //     double pot_act_tmp[N_MAX];
+// // // // // // //     double3 a_act_tmp[N_MAX];
+// // // // // // //     double3 adot_act_tmp[N_MAX];
+
     g6_set_ti(clusterid, t);
 
     int ni = n_act; // TODO why is this needed?
@@ -1067,12 +1016,12 @@ void calc_self_grav(double t, double eps2, double &g6_calls)
 
     /* Store the new value of the local partial force etc... */
 
-    for (int i=0; i<n_act; i++) {
+// // // // // // // //     for (int i=0; i<n_act; i++) {
-        int iii = ind_act[i];
+// // // // // // // //         int iii = ind_act[i];
-        pot_act_tmp_loc[iii] = pot_act_tmp[i];
+// // // // // // // //         pot_act_tmp_loc[iii] = pot_act_tmp[i];
-        a_act_tmp_loc[iii] = a_act_tmp[i];
+// // // // // // // //         a_act_tmp_loc[iii] = a_act_tmp[i];
-        adot_act_tmp_loc[iii] = adot_act_tmp[i];
+// // // // // // // //         adot_act_tmp_loc[iii] = adot_act_tmp[i];
-    } /* i */
+// // // // // // // //     } /* i */
 
 #ifdef TIMING
     get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
@@ -1337,7 +1286,7 @@ DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
 #endif        // EXTPOT
 }
 
-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) // TODO how do we know N?
+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[])
 {
     // TODO maybe use static variables here for the previous step energy?
     double E_pot = 0.0;
@@ -1380,7 +1329,7 @@ void energy_contr(const double time_cur, const double timesteps, const double n_
     rcm_sum += rcm_mod;
     vcm_sum += vcm_mod;
 
-    mom = {0, 0, 0};
+    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] );
@@ -1405,7 +1354,7 @@ void energy_contr(const double time_cur, const double timesteps, const double n_
     } else {
         /* read the E_tot_0 + etc... */
         if (skip_con == 0) {
-            inp = fopen("contr.con","r");
+            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,
@@ -1428,7 +1377,7 @@ void energy_contr(const double time_cur, const double timesteps, const double n_
 
     fflush(stdout);
 
-    out = fopen("contr.dat", "a");
+    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,
@@ -1438,7 +1387,7 @@ void energy_contr(const double time_cur, const double timesteps, const double n_
             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");
+    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,
@@ -1474,6 +1423,45 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
     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;
+
 
     /* INIT the rand() !!! */
     srand(19640916);                 /* it is just my birthday :-) */
@@ -1571,7 +1559,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
 
         /* read the global data for particles to the rootRank */
 
-        read_data(&diskstep, &N, &time_cur, ind, m, x, v);
+        read_data(inp_fname, &diskstep, &N, &time_cur, ind, m, x, v);
 
         /* possible coordinate & velocity limits for ALL particles !!! */
 
@@ -1851,10 +1839,14 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
 
     } /* i */
 
-    // NOTE this is where calc_self_grav_zero() used to be. Hopefully the copying from *_act to *_act_new is just a temporary measure.
+    // NOTE this is where calc_self_grav_zero() used to be.
-    std::copy(&x_act[0][0], &x_act[0][0]+3*N, &x_act_new[0][0]);
+    calc_self_grav(time_cur, eps2, g6_calls, n_loc,
-    std::copy(&v_act[0][0], &v_act[0][0]+3*N, &v_act_new[0][0]);
+                   n_act, ind_act,
-    calc_self_grav(time_cur, eps2, g6_calls);
+                   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);
@@ -1965,7 +1957,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
 
     /* 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);
+        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);
     } /* if (myRank == rootRank) */
 
 #ifdef ETICS_DUMP
@@ -2062,10 +2054,10 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
     if (myRank == rootRank) {
 
         /* Write BH data... */
-        if (config->live_smbh_output) write_bh_data(time_cur, m, x, v);
+        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, m, x, v);
+        if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur, N, m, x, v);
 
     } /* if (myRank == rootRank) */
 
@@ -2261,7 +2253,13 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
         DT_ACT_PRED += (CPU_tmp_user - CPU_tmp_user0);
 #endif
 
-        calc_self_grav(min_t, eps2, g6_calls);
+        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 */
 
@@ -2580,10 +2578,10 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
         if (time_cur >= t_bh) {
             if (myRank == rootRank) {
                 /* Write BH data... */
-                if (config->live_smbh_output) write_bh_data(time_cur, m, x, v);
+                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, m, x, v);
+                if (config->live_smbh_neighbor_output) write_bh_nb_data(time_cur, N, m, x, v);
 
             } /* if (myRank == rootRank) */
 
@@ -2593,7 +2591,7 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0,
         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);
+                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);
 
                 /* write cont data */