diff --git a/phigrape.conf b/phigrape.conf index 337b154..2414d1e 100644 --- a/phigrape.conf +++ b/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_b_bulge = 20 -ext_m_disk = 1e6 -ext_a_disk = 80 -ext_b_disk = 70 -ext_m_halo_plummer = 1e6 -ext_b_halo_plummer = 200. +ext_m_bulge = 0 +ext_b_bulge = 0 +ext_m_disk = 1e9 +ext_a_disk = 0.8 +ext_b_disk = 0.2 +ext_m_halo_plummer = 0 +ext_b_halo_plummer = 0 #################################### diff --git a/phigrape.cpp b/phigrape.cpp index f46dc36..e4153a1 100644 --- a/phigrape.cpp +++ b/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 calc_plummer(const int n_act, const double3 x[], const double3 v[], double pot[], double3 a[], double3 adot[], int mode) + void apply(const int n_act, const double3 x[], const double3 v[], double pot[], double3 a[], double3 adot[]) { - /* 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; iset_coordinates(x[i], v[i]); + this->calc_gravity(); + pot[i] += potential; + a[i] += acceleration; + adot[i] += jerk; } } - 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; - double z2_tmp = x[i][2]*x[i][2] + b_disk*b_disk; - 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; - + this->x = x; + this->v = v; } -//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_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) - 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 + std::fill(pot_act_ext, pot_act_ext+n_act, 0.); - - -for (int i=0; iis_active) + component->apply(n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new); } #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; iext_m_bulge/config->norm_mass; - external_gravity.b_bulge = config->ext_b_bulge/config->norm_length*1000; - external_gravity.m_disk = config->ext_m_disk/config->norm_mass; - external_gravity.a_disk = config->ext_a_disk/config->norm_length*1000; - external_gravity.b_disk = config->ext_b_disk/config->norm_length*1000; - external_gravity.m_halo = config->ext_m_halo_plummer/config->norm_mass; - external_gravity.b_halo = config->ext_b_halo_plummer/config->norm_length*1000; + double normalization_mass=1, normalization_length=1; + if (config->ext_units_physical) { + normalization_mass = 1/config->norm_mass; + normalization_length = 1000/config->norm_length; + } + 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); - if (external_gravity.m_bulge > 0) - printf("m_bulge = %.4E b_bulge = %.4E\n", external_gravity.m_bulge, external_gravity.b_bulge); - if (external_gravity.m_disk > 0) - printf("m_disk = %.4E a_disk = %.4E b_disk = %.4E\n", external_gravity.m_disk, external_gravity.a_disk, external_gravity.b_disk); - if (external_gravity.m_halo > 0) - printf("m_halo = %.4E b_halo = %.4E\n", external_gravity.m_halo, external_gravity.b_halo); - + std::vector external_gravity_components; + external_gravity_components.push_back(&ext_bulge); + external_gravity_components.push_back(&ext_disk); + external_gravity_components.push_back(&ext_halo_plummer); + + 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_zero();*/*/*/*/*/*/ - calc_ext_grav(n_act, x_act, v_act, a, adot); -#endif + calc_ext_grav(external_gravity_components, n_act, x_act, v_act, pot_ext, a, adot); /* 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(n_act, x_act_new, v_act_new, a_act_new, adot_act_new); -#endif + calc_ext_grav(external_gravity_components, n_act, x_act_new, v_act_new, pot_act_ext, a_act_new, adot_act_new); /* 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];