New approach to external potential

phigrape.conf

@@ -40,13 +40,13 @@ norm_length = 10
 
 # Notice that if physical units are used, the "a" and "b" parameters are in kiloparsec, not parsec!
 
-ext_m_bulge = 1e5
+ext_m_bulge = 0
-ext_b_bulge = 20
+ext_b_bulge = 0
-ext_m_disk = 1e6
+ext_m_disk = 1e9
-ext_a_disk = 80
+ext_a_disk = 0.8
-ext_b_disk = 70
+ext_b_disk = 0.2
-ext_m_halo_plummer = 1e6
+ext_m_halo_plummer = 0
-ext_b_halo_plummer = 200.
+ext_b_halo_plummer = 0
 
 
 ####################################

phigrape.cpp

@@ -58,7 +58,6 @@ Config *config;
 
 #define NORM            // Physical normalization
 
-#define EXTPOT            // external potential (BH? or galactic?)
 //#define EXTPOT_BH        // BH - usually NB units
 
 
@@ -307,7 +306,6 @@ double eps_BH=0.0;
 
 /* external potential... */
 
-#ifdef EXTPOT
 
 #ifdef EXTPOT_GAL_DEH
 double m_ext, r_ext, g_ext,
@@ -329,7 +327,6 @@ double  r2, rv_ij,
 
 double pot_ext[N_MAX], pot_act_ext[N_MAX];    // for all EXTPOT
 
-#endif
 
 int    i_bh, i_bh1, i_bh2,
     num_bh = 0, num_bh1 = 0, num_bh2 = 0;
@@ -625,70 +622,85 @@ void write_bh_nb_data(double time_cur, int N, double m[], double3 x[], double3 v
 
 class External_gravity {
 public:
-    External_gravity() {}
+    void apply(const int n_act, const double3 x[], const double3 v[], double pot[], double3 a[], double3 adot[])
-    void calc_plummer(const int n_act, const double3 x[], const double3 v[], double pot[], double3 a[], double3 adot[], int mode)
     {
-        /* This funcion calculates gravity due to a smeared point source using Plummer softening. The `mode` parameter determines whether we use the BH, bulge, or halo. */ // TODO change mode to class enum?
-
-        double m, b;
-        switch (mode) {
-            case 0: m=m_bh;    b=b_bh;    break;
-            case 1: m=m_bulge; b=b_bulge; break;
-            case 2: m=m_halo;  b=b_halo;  break;
-            default:
-                throw std::runtime_error("Unknown mode in external_gravity.calc_plummer");
-        }
         for (int i=0; i<n_act; i++) {
-            double r2 = SQR(b);
+            this->set_coordinates(x[i], v[i]);
-            r2 += x[i].norm2();
+            this->calc_gravity();
-            double r = sqrt(r2);
+            pot[i]  += potential;
-            double rv_ij = v[i]*x[i];
+            a[i]    += acceleration;
-            double tmp = m / r;
+            adot[i] += jerk;
-            pot[i] -= tmp;
-            tmp /= r2;
-            a[i] -= tmp * x[i];
-            adot[i] -= tmp * (v[i] - 3*rv_ij * x[i]/r2);
         }
     }
-    void calc_disk(const int n_act, double3 x[], double3 v[], double pot[], double3 a[], double3 adot[]) //TODO change x and v to const
+    virtual void calc_gravity() = 0;
+    bool is_active;
+protected:
+    double potential;
+    double3 acceleration, jerk;
+    double3 x, v;
+    void set_coordinates(double3 x, double3 v)
     {
-        int i = 0;
+        this->x = x;
-        double z2_tmp = x[i][2]*x[i][2] + b_disk*b_disk;
+        this->v = v;
-        double z_tmp  = sqrt(z2_tmp);
-        double r2_tmp = x[i][0]*x[i][0] + x[i][1]*x[i][1] + SQR(z_tmp + a_disk);
-        double r_tmp  = sqrt(r2_tmp);
-        pot[i] -= m_disk / r_tmp;
-        double tmp = m_disk / (r2_tmp*r_tmp);
-        a[i][0] -= tmp * x_ij;
-        a[i][1] -= tmp * y_ij;
-        a[i][2] -= tmp * z_ij * (z_tmp + a_disk)/z_tmp;
-        tmp = m_disk / (z_tmp*r2_tmp*r2_tmp*r_tmp);
-        adot[i][0] += tmp * (- v[i][0]*z_tmp*r2_tmp
-                      + 3*x[i][0]*v[i][0]*x[i][0]*z_tmp
-                      + 3*x[i][0]*v[i][1]*x[i][1]*z_tmp
-                      + 3*x[i][0]*v[i][2]*x[i][2]*SQR(z_tmp + a_disk));
-        adot[i][1] += tmp * (- v[i][1]*z_tmp*r2_tmp
-                      + 3*x[i][1]*v[i][0]*x[i][0]*z_tmp
-                      + 3*x[i][1]*v[i][1]*x[i][1]*z_tmp
-                      + 3*x[i][1]*v[i][2]*x[i][2]*SQR(z_tmp + a_disk));
-        adot[i][2] += tmp * (- v[i][2]*(z_tmp + a_disk)*(x[i][0]*(z2_tmp*z_tmp + a_disk*SQR(b_disk)) +
-                      x[i][1]*(z2_tmp*z_tmp + a_disk*SQR(b_disk)) -
-                      (2*a_disk*(SQR(x[i][2]) - SQR(b_disk))*z_tmp +
-                      2*SQR(x[i][2]*x[i][2]) +
-                      SQR(b_disk)*SQR(x[i][2]) -
-                      SQR(b_disk)*( SQR(a_disk) + SQR(b_disk))))
-                      + 3*v[i][0]*x[i][0]*x[i][2]*z2_tmp*(z_tmp + a_disk)
-                      + 3*v[i][1]*x[i][1]*x[i][2]*z2_tmp*(z_tmp + a_disk)) / z2_tmp;
-
     }
-//private: //TODO make private
-    double m_bh,            b_bh,
-           m_bulge,         b_bulge,
-           m_disk,  a_disk, b_disk,
-           m_halo,          b_halo;
 };
 
-External_gravity external_gravity;
+class Plummer : public External_gravity {
+public:
+    Plummer(double m, double b) : m(m), b(b) {is_active=(m>0);}
+private:
+    void calc_gravity()
+    {
+        double r2 = SQR(b);
+        r2 += x.norm2();
+        double r = sqrt(r2);
+        double rv_ij = v*x;
+        double tmp = m / r;
+        potential = -tmp;
+        tmp /= r2;
+        acceleration = - tmp * x;
+        jerk = - tmp * (v - 3*rv_ij * x/r2);
+    }
+    double m, b;
+};
+
+class Miyamoto_Nagai : public External_gravity {
+public:
+    Miyamoto_Nagai(double m, double a, double b) : m(m), a(a), b(b) {is_active=(m>0);}
+private:
+    void calc_gravity()
+    {
+        double  x_ij=x[0],  y_ij=x[1],  z_ij=x[2];
+        double vx_ij=v[0], vy_ij=v[1], vz_ij=v[2];
+        auto z2_tmp = z_ij*z_ij + SQR(b);
+        auto z_tmp  = sqrt(z2_tmp);
+        auto r2_tmp = x_ij*x_ij + y_ij*y_ij + SQR(z_tmp + a);
+        auto r_tmp  = sqrt(r2_tmp);
+        potential = - m / r_tmp;
+        auto tmp = m / (r2_tmp*r_tmp);
+        acceleration[0] = - tmp * x_ij;
+        acceleration[1] = - tmp * y_ij;
+        acceleration[2] = - tmp * z_ij * (z_tmp + a)/z_tmp;
+        tmp = m / (z_tmp*r2_tmp*r2_tmp*r_tmp);
+        jerk[0] = tmp * (- vx_ij*z_tmp*r2_tmp
+                + 3.0*x_ij*vx_ij*x_ij*z_tmp
+                + 3.0*x_ij*vy_ij*y_ij*z_tmp
+                + 3.0*x_ij*vz_ij*z_ij*SQR(z_tmp + a));
+        jerk[1] = tmp * (- vy_ij*z_tmp*r2_tmp
+                + 3.0*y_ij*vx_ij*x_ij*z_tmp
+                + 3.0*y_ij*vy_ij*y_ij*z_tmp
+                + 3.0*y_ij*vz_ij*z_ij*SQR(z_tmp + a));
+        jerk[2] = tmp * (- vz_ij*(z_tmp + a)*(x_ij*x_ij*(z2_tmp*z_tmp + a*SQR(b)) +
+                y_ij*y_ij*(z2_tmp*z_tmp + a*SQR(b)) -
+                (2.0*a*(SQR(z_ij) - SQR(b))*z_tmp +
+                2.0*SQR(z_ij*z_ij) +
+                SQR(b)*SQR(z_ij) -
+                SQR(b)*(SQR(a) + SQR(b))))
+                + 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a)
+                + 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a)) / z2_tmp;
+    }
+    double m, a, b;
+};
 
 void calc_self_grav(double t, double eps2, double &g6_calls, int n_loc,
                     int n_act, int ind_act[],
@@ -728,27 +740,18 @@ void calc_self_grav(double t, double eps2, double &g6_calls, int n_loc,
 #endif
 }
 
-void calc_ext_grav(int n_act, double3 *x_act_new, double3 *v_act_new, double3 *a_act_new, double3* adot_act_new)
+void calc_ext_grav(std::vector<External_gravity*> &external_gravity_components, int n_act, double3 *x_act_new, double3 *v_act_new, double *pot_act_ext, double3 *a_act_new, double3* adot_act_new)
 {
 
-#ifdef EXTPOT
-
 /* Define the external potential for all active particles on all nodes */
 
-int ni = n_act;
+int ni = n_act; // TODO redundant?
 
-    if (external_gravity.m_bulge > 0)
+    std::fill(pot_act_ext, pot_act_ext+n_act, 0.);
-        external_gravity.calc_plummer(ni, x_act_new,v_act_new, pot_act_ext, a_act_new, adot_act_new, 1); // Bulge
-    if (external_gravity.m_disk > 0)
-        external_gravity.calc_disk(ni, x_act_new,v_act_new, pot_act_ext, a_act_new, adot_act_new); // Disk
-    if (external_gravity.m_halo > 0)
-        external_gravity.calc_plummer(ni, x_act_new,v_act_new, pot_act_ext, a_act_new, adot_act_new, 2); // Halo
 
-
+    for (auto component : external_gravity_components) {
-
+        if (component->is_active)
-for (int i=0; i<ni; i++)
+            component->apply(n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new);
-    {
-    pot_act_ext[i] = 0.0;
     }
 
 #ifdef TIMING
@@ -880,7 +883,6 @@ get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
 DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
 #endif
 
-#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, int N, double m[], double3 x[], double3 v[], double pot[])
@@ -895,9 +897,8 @@ void energy_contr(const double time_cur, const double timesteps, const double n_
     E_kin *= 0.5;
 
     double E_pot_ext = 0.0;
-#ifdef EXTPOT
+
     for (int i=0; i<N; i++) E_pot_ext += m[i]*pot_ext[i];
-#endif
 
     if (timesteps == 0.0) { //TODO what is this thing?
         rcm_sum = 0.0;
@@ -1108,8 +1109,6 @@ int main(int argc, char *argv[])
         */
 
         /* Read the data for EXTernal POTential */
-#ifdef EXTPOT
-//         auto inp = fopen("phi-GRAPE.ext","r");
 
 #ifdef EXTPOT_BH
         fscanf(inp,"%lE %lE", &m_bh, &b_bh);
@@ -1127,9 +1126,6 @@ int main(int argc, char *argv[])
         r2_halo = SQR(r_halo);
 #endif
 
-//         fclose(inp);
-#endif        // EXTPOT
-
         /* read the global data for particles to the rootRank */
 
         read_data(inp_fname, &diskstep, &N, &time_cur, ind, m, x, v);
@@ -1154,7 +1150,6 @@ int main(int argc, char *argv[])
         printf("\n");
 #endif
 
-#ifdef EXTPOT
         printf("External Potential: \n");
         printf("\n");
 
@@ -1170,10 +1165,6 @@ int main(int argc, char *argv[])
         printf("v_halo = %.6E r_halo = %.6E \n", v_halo, r_halo);
 #endif
 
-        printf("\n");
-#endif
-
-        fflush(stdout);
 
         if ((diskstep == 0) && (time_cur == 0.0)) {
             out = fopen("contr.dat", "w");
@@ -1225,8 +1216,6 @@ int main(int argc, char *argv[])
     MPI_Bcast(&t_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
 #endif        // NORM
 
-#ifdef EXTPOT
-
 #ifdef EXTPOT_BH
     MPI_Bcast(&m_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&b_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
@@ -1246,26 +1235,28 @@ int main(int argc, char *argv[])
     MPI_Bcast(&r2_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
 #endif
 
-#endif        // EXTPOT
-
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
-    external_gravity.m_bulge = config->ext_m_bulge/config->norm_mass;
+    double normalization_mass=1, normalization_length=1;
-    external_gravity.b_bulge = config->ext_b_bulge/config->norm_length*1000;
+    if (config->ext_units_physical) {
-    external_gravity.m_disk  = config->ext_m_disk/config->norm_mass;
+        normalization_mass = 1/config->norm_mass;
-    external_gravity.a_disk  = config->ext_a_disk/config->norm_length*1000;
+        normalization_length = 1000/config->norm_length;
-    external_gravity.b_disk  = config->ext_b_disk/config->norm_length*1000;
+    }
-    external_gravity.m_halo  = config->ext_m_halo_plummer/config->norm_mass;
+    Plummer ext_bulge(config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length);
-    external_gravity.b_halo  = config->ext_b_halo_plummer/config->norm_length*1000;
+    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);
 
-    if (external_gravity.m_bulge > 0)
+    std::vector<External_gravity*> external_gravity_components;
-        printf("m_bulge = %.4E                        b_bulge = %.4E\n", external_gravity.m_bulge, external_gravity.b_bulge);
+    external_gravity_components.push_back(&ext_bulge);
-    if (external_gravity.m_disk > 0)
+    external_gravity_components.push_back(&ext_disk);
-        printf("m_disk  = %.4E  a_disk  = %.4E  b_disk  = %.4E\n", external_gravity.m_disk, external_gravity.a_disk, external_gravity.b_disk);
+    external_gravity_components.push_back(&ext_halo_plummer);
-    if (external_gravity.m_halo > 0)
+
-        printf("m_halo  = %.4E                        b_halo  = %.4E\n", external_gravity.m_halo, external_gravity.b_halo);
+    if (config->ext_m_bulge > 0)        printf("m_bulge = %.4E                        b_bulge = %.4E\n", config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length);
-    
+    if (config->ext_m_disk > 0)         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 (config->ext_m_halo_plummer > 0) printf("m_halo  = %.4E                        b_halo = %.4E\n", config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length);
+    printf("\n");
+    fflush(stdout);
 
     eta_s  = eta/ETA_S_CORR;
     eta_bh = eta/ETA_BH_CORR;
@@ -1502,10 +1493,7 @@ int main(int argc, char *argv[])
         }
     }
 
-#ifdef EXTPOT
+    calc_ext_grav(external_gravity_components, n_act, x_act, v_act, pot_ext, a, adot);
-//     /*/*/*/*/*/*calc_ext_grav_zero();*/*/*/*/*/*/
-    calc_ext_grav(n_act, x_act, v_act, a, adot);
-#endif
 
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
@@ -1777,9 +1765,7 @@ int main(int argc, char *argv[])
 
             pot_act[i] = pot[iii];
 
-#ifdef EXTPOT
             pot_act_ext[i] = pot_ext[iii];
-#endif
 
             a_act[i] = a[iii];
             adot_act[i] = adot[iii];
@@ -1910,9 +1896,7 @@ int main(int argc, char *argv[])
             }
         }
 
-#ifdef EXTPOT
+        calc_ext_grav(external_gravity_components, n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new);
-        calc_ext_grav(n_act, x_act_new, v_act_new, a_act_new, adot_act_new);
-#endif
 
         /* correct the active particles positions etc... on all the nodes */
 
@@ -2080,9 +2064,7 @@ int main(int argc, char *argv[])
 
             pot[iii] = pot_act[i];
 
-#ifdef EXTPOT
             pot_ext[iii] = pot_act_ext[i];
-#endif
 
             a[iii] = a_act[i];
             adot[iii] = adot_act[i];