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

Added logarithmic external potential

Browse source
This commit is contained in:
Yohai Meiron 2020-03-25 21:42:29 -04:00
parent 8cba5b5e1d
commit 131642aa71
5 changed files with 51 additions and 84 deletions
Download patch file Download diff file Expand all files Collapse all files

100
phigrape.cpp
View file

@ -62,7 +62,6 @@ Config *config;
//#define EXTPOT_GAL_DEH // Dehnen Galactic - usually physical units
//#define EXTPOT_GAL_LOG // Log Galactic - usually physical units
#ifdef ETICS
#include "grapite.h"
@ -312,12 +311,6 @@ double m_ext, r_ext, g_ext,
tmp_r2, tmp_r3, dum, dum2, dum3, dum_g, tmp_g;
#endif
#ifdef EXTPOT_GAL_LOG
double v_halo, r_halo,
v2_halo, r2_halo, r2_r2_halo,
x2_ij, y2_ij, z2_ij;
#endif
#ifdef EXTPOT_BH
double r2, rv_ij,
m_bh, b_bh, eps_bh; // NOTE different from eps_BH
@ -591,6 +584,24 @@ private:
double m, a, b;
};
class Logarithmic_halo : public External_gravity {
public:
Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0);}
private:
void calc_gravity()
{
auto r2 = x.norm2();
auto rv_ij = 2.0*(v*x);
auto r2_r2_halo = (r2 + r2_halo);
potential = - 0.5*v2_halo * log(1.0 + r2/r2_halo);
auto tmp = v2_halo/r2_r2_halo;
acceleration = - tmp * x;
tmp /= (r2 + r2_halo);
jerk = tmp * (rv_ij * x - r2_r2_halo * v);
}
double v2_halo, r2_halo;
};
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[],
@ -688,44 +699,6 @@ for (i=0; i<ni; i++)
} /* i */
#endif
/* For simple Log potential... */
#ifdef EXTPOT_GAL_LOG
for (i=0; i<ni; i++)
{
x_ij = x_act_new[i][0];
y_ij = x_act_new[i][1];
z_ij = x_act_new[i][2];
vx_ij = v_act_new[i][0];
vy_ij = v_act_new[i][1];
vz_ij = v_act_new[i][2];
r2 = SQR(x_ij) + SQR(y_ij) + SQR(z_ij);
rv_ij = 2.0 * ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
r2_r2_halo = (r2 + r2_halo);
pot_act_ext[i] -= 0.5*v2_halo * log(1.0 + r2/r2_halo);
tmp = v2_halo/r2_r2_halo;
a_act_new[i][0] -= tmp * x_ij;
a_act_new[i][1] -= tmp * y_ij;
a_act_new[i][2] -= tmp * z_ij;
tmp /= (r2 + r2_halo);
adot_act_new[i][0] += tmp * (rv_ij * x_ij - r2_r2_halo * vx_ij);
adot_act_new[i][1] += tmp * (rv_ij * y_ij - r2_r2_halo * vy_ij);
adot_act_new[i][2] += tmp * (rv_ij * z_ij - r2_r2_halo * vz_ij);
} /* i */
#endif
/* For simple Plummer potential... */
@ -1022,14 +995,6 @@ int main(int argc, char *argv[])
fscanf(inp,"%lE %lE %lE", &m_ext, &r_ext, &g_ext);
#endif
#ifdef EXTPOT_GAL_LOG
fscanf(inp,"%lE %lE", &v_halo, &r_halo);
v_halo *= (km/v_norm); r_halo *= (pc/r_norm);
v2_halo = SQR(v_halo);
r2_halo = SQR(r_halo);
#endif
/* read the global data for particles to the rootRank */
read_data(inp_fname, &diskstep, &N, &time_cur, ind, m, x, v);
@ -1065,11 +1030,6 @@ int main(int argc, char *argv[])
printf("m_ext = %.6E r_ext = %.6E \t g_ext = %.3E \n", m_ext, r_ext, g_ext);
#endif
#ifdef EXTPOT_GAL_LOG
printf("v_halo = %.6E r_halo = %.6E \n", v_halo, r_halo);
#endif
if ((diskstep == 0) && (time_cur == 0.0)) {
out = fopen("contr.dat", "w");
fclose(out);
@ -1131,34 +1091,30 @@ int main(int argc, char *argv[])
MPI_Bcast(&g_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
#endif
#ifdef EXTPOT_GAL_LOG
MPI_Bcast(&v_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&r_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&v2_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
MPI_Bcast(&r2_halo, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
#endif
/* Wait to all processors to finish his works... */
MPI_Barrier(MPI_COMM_WORLD);
double normalization_mass=1, normalization_length=1;
double normalization_mass=1, normalization_length=1, normalization_velocity=1;
if (config->ext_units_physical) {
normalization_mass = 1/config->norm_mass;
normalization_length = 1000/config->norm_length;
normalization_mass = 1/config->unit_mass;
normalization_length = 1000/config->unit_length;
normalization_velocity = 1.52484071426404437233e+01*sqrt(config->unit_length/config->unit_mass);
}
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);
Logarithmic_halo ext_log_halo(config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length);
std::vector<External_gravity*> 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);
external_gravity_components.push_back(&ext_log_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);
if (ext_bulge.is_active) printf("m_bulge = %.4E b_bulge = %.4E\n", config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length);
if (ext_disk.is_active) 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 (ext_halo_plummer.is_active) printf("m_halo = %.4E b_halo = %.4E\n", config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length);
if (ext_log_halo.is_active) printf("v_halo = %.6E r_halo = %.6E \n", config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length);
printf("\n");
fflush(stdout);