Skip to content

GitLab

Commit df794475 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'gravitational_potential' into 'master' 

Gravitational potential

See merge request !512
parents ca44756d 215d887d
Hide whitespace changes
Inline Side-by-side
.gitignore
View file @ df794475
......@@ -105,6 +105,8 @@ tests/testDump
tests/testLogger
tests/benchmarkInteractions
tests/testGravityDerivatives
tests/testPotentialSelf
tests/testPotentialPair
theory/latex/swift.pdf
theory/SPH/Kernels/kernels.pdf
......
examples/plot_gravity_checks.py
View file @ df794475
......@@ -13,7 +13,7 @@ params = {'axes.labelsize': 14,
'xtick.labelsize': 14,
'ytick.labelsize': 14,
'text.usetex': True,
'figure.figsize': (10, 10),
'figure.figsize': (12, 10),
'figure.subplot.left' : 0.06,
'figure.subplot.right' : 0.99 ,
'figure.subplot.bottom' : 0.06 ,
......@@ -60,11 +60,13 @@ data = np.loadtxt('gravity_checks_exact_step%d.dat'%step)
exact_ids = data[:,0]
exact_pos = data[:,1:4]
exact_a = data[:,4:7]
exact_pot = data[:,7]
# Sort stuff
sort_index = np.argsort(exact_ids)
exact_ids = exact_ids[sort_index]
exact_pos = exact_pos[sort_index, :]
exact_a = exact_a[sort_index, :]
exact_pot = exact_pot[sort_index]
exact_a_norm = np.sqrt(exact_a[:,0]**2 + exact_a[:,1]**2 + exact_a[:,2]**2)
# Start the plot
......@@ -166,12 +168,14 @@ for i in range(num_order):
ids = data[:,0]
pos = data[:,1:4]
a_grav = data[:, 4:7]
pot = data[:, 7]
# Sort stuff
sort_index = np.argsort(ids)
ids = ids[sort_index]
pos = pos[sort_index, :]
a_grav = a_grav[sort_index, :]
pot = pot[sort_index]
# Cross-checks
if not np.array_equal(exact_ids, ids):
......@@ -185,6 +189,7 @@ for i in range(num_order):
# Compute the error norm
diff = exact_a - a_grav
diff_pot = exact_pot - pot
norm_diff = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)
......@@ -192,73 +197,87 @@ for i in range(num_order):
error_x = abs(diff[:,0]) / exact_a_norm
error_y = abs(diff[:,1]) / exact_a_norm
error_z = abs(diff[:,2]) / exact_a_norm
error_pot = abs(diff_pot) / abs(exact_pot)
# Bin the error
norm_error_hist,_ = np.histogram(norm_error, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
error_x_hist,_ = np.histogram(error_x, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
error_y_hist,_ = np.histogram(error_y, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
error_z_hist,_ = np.histogram(error_z, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
error_pot_hist,_ = np.histogram(error_pot, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
norm_median = np.median(norm_error)
median_x = np.median(error_x)
median_y = np.median(error_y)
median_z = np.median(error_z)
median_pot = np.median(error_pot)
norm_per99 = np.percentile(norm_error,99)
per99_x = np.percentile(error_x,99)
per99_y = np.percentile(error_y,99)
per99_z = np.percentile(error_z,99)
per99_pot = np.percentile(error_pot, 99)
norm_max = np.max(norm_error)
max_x = np.max(error_x)
max_y = np.max(error_y)
max_z = np.max(error_z)
max_pot = np.max(error_pot)
print "Order %d ---- "%order[i]
print "Norm: median= %f 99%%= %f max= %f"%(norm_median, norm_per99, norm_max)
print "X : median= %f 99%%= %f max= %f"%(median_x, per99_x, max_x)
print "Y : median= %f 99%%= %f max= %f"%(median_y, per99_y, max_y)
print "Z : median= %f 99%%= %f max= %f"%(median_z, per99_z, max_z)
print "Pot : median= %f 99%%= %f max= %f"%(median_pot, per99_pot, max_pot)
print ""
plt.subplot(221)
plt.semilogx(bins, norm_error_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(norm_median, norm_per99), ha="left", va="top", color=cols[i])
plt.subplot(222)
plt.subplot(231)
plt.semilogx(bins, error_x_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_x, per99_x), ha="left", va="top", color=cols[i])
plt.subplot(223)
plt.subplot(232)
plt.semilogx(bins, error_y_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_y, per99_y), ha="left", va="top", color=cols[i])
plt.subplot(224)
plt.subplot(233)
plt.semilogx(bins, error_z_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_z, per99_z), ha="left", va="top", color=cols[i])
plt.subplot(234)
plt.semilogx(bins, norm_error_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(norm_median, norm_per99), ha="left", va="top", color=cols[i])
plt.subplot(235)
plt.semilogx(bins, error_pot_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_pot, per99_pot), ha="left", va="top", color=cols[i])
count += 1
plt.subplot(221)
plt.xlabel("$|\delta \overrightarrow{a}|/|\overrightarrow{a}_{exact}|$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,2.5)
plt.legend(loc="upper left")
plt.subplot(222)
plt.subplot(231)
plt.xlabel("$\delta a_x/|\overrightarrow{a}_{exact}|$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
#plt.legend(loc="center left")
plt.subplot(223)
plt.subplot(232)
plt.xlabel("$\delta a_y/|\overrightarrow{a}_{exact}|$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
#plt.legend(loc="center left")
plt.subplot(224)
plt.subplot(233)
plt.xlabel("$\delta a_z/|\overrightarrow{a}_{exact}|$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
plt.subplot(234)
plt.xlabel("$|\delta \overrightarrow{a}|/|\overrightarrow{a}_{exact}|$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,2.5)
plt.legend(loc="upper left")
plt.subplot(235)
plt.xlabel("$\delta \phi/\phi_{exact}$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
#plt.legend(loc="center left")
......
src/gravity.c
View file @ df794475
......@@ -54,6 +54,7 @@ struct exact_force_data {
static float fewald_x[Newald + 1][Newald + 1][Newald + 1];
static float fewald_y[Newald + 1][Newald + 1][Newald + 1];
static float fewald_z[Newald + 1][Newald + 1][Newald + 1];
static float potewald[Newald + 1][Newald + 1][Newald + 1];
/* Factor used to normalize the access to the Ewald table */
float ewald_fac;
......@@ -83,6 +84,8 @@ void gravity_exact_force_ewald_init(double boxSize) {
const float factor_exp1 = 2.f * alpha / sqrt(M_PI);
const float factor_exp2 = -M_PI * M_PI / alpha2;
const float factor_sin = 2.f * M_PI;
const float factor_cos = 2.f * M_PI;
const float factor_pot = M_PI / alpha2;
const float boxSize_inv2 = 1.f / (boxSize * boxSize);
/* Ewald factor to access the table */
......@@ -92,6 +95,9 @@ void gravity_exact_force_ewald_init(double boxSize) {
bzero(fewald_x, (Newald + 1) * (Newald + 1) * (Newald + 1) * sizeof(float));
bzero(fewald_y, (Newald + 1) * (Newald + 1) * (Newald + 1) * sizeof(float));
bzero(fewald_z, (Newald + 1) * (Newald + 1) * (Newald + 1) * sizeof(float));
bzero(potewald, (Newald + 1) * (Newald + 1) * (Newald + 1) * sizeof(float));
potewald[0][0][0] = 2.8372975f;
/* Compute the values in one of the octants */
for (int i = 0; i <= Newald; ++i) {
......@@ -114,6 +120,7 @@ void gravity_exact_force_ewald_init(double boxSize) {
float f_x = r_x * r_inv3;
float f_y = r_y * r_inv3;
float f_z = r_z * r_inv3;
float pot = r_inv + factor_pot;
for (int n_i = -4; n_i <= 4; ++n_i) {
for (int n_j = -4; n_j <= 4; ++n_j) {
......@@ -130,15 +137,17 @@ void gravity_exact_force_ewald_init(double boxSize) {
const float r_tilde_inv3 =
r_tilde_inv * r_tilde_inv * r_tilde_inv;
const float val =
erfcf(alpha * r_tilde) +
factor_exp1 * r_tilde * expf(-alpha2 * r_tilde2);
const float val_pot = erfcf(alpha * r_tilde);
const float val_f =
val_pot + factor_exp1 * r_tilde * expf(-alpha2 * r_tilde2);
/* First correction term */
const float f = val * r_tilde_inv3;
const float f = val_f * r_tilde_inv3;
f_x -= f * d_x;
f_y -= f * d_y;
f_z -= f * d_z;
pot -= val_pot * r_tilde_inv;
}
}
}
......@@ -149,16 +158,23 @@ void gravity_exact_force_ewald_init(double boxSize) {
const float h2 = h_i * h_i + h_j * h_j + h_k * h_k;
if (h2 == 0.f) continue;
const float h2_inv = 1.f / (h2 + FLT_MIN);
const float h_dot_x = h_i * r_x + h_j * r_y + h_k * r_z;
const float val = 2.f * h2_inv * expf(h2 * factor_exp2) *
sinf(factor_sin * h_dot_x);
const float common = h2_inv * expf(h2 * factor_exp2);
const float val_pot =
(float)M_1_PI * common * cosf(factor_cos * h_dot_x);
const float val_f = 2.f * common * sinf(factor_sin * h_dot_x);
/* Second correction term */
f_x -= val * h_i;
f_y -= val * h_j;
f_z -= val * h_k;
f_x -= val_f * h_i;
f_y -= val_f * h_j;
f_z -= val_f * h_k;
pot -= val_pot;
}
}
}
......@@ -167,6 +183,7 @@ void gravity_exact_force_ewald_init(double boxSize) {
fewald_x[i][j][k] = f_x;
fewald_y[i][j][k] = f_y;
fewald_z[i][j][k] = f_z;
potewald[i][j][k] = pot;
}
}
}
......@@ -196,6 +213,11 @@ void gravity_exact_force_ewald_init(double boxSize) {
H5Dwrite(h_data, H5T_NATIVE_FLOAT, h_space, H5S_ALL, H5P_DEFAULT,
&(fewald_z[0][0][0]));
H5Dclose(h_data);
h_data = H5Dcreate(h_file, "Ewald_pot", H5T_NATIVE_FLOAT, h_space,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
H5Dwrite(h_data, H5T_NATIVE_FLOAT, h_space, H5S_ALL, H5P_DEFAULT,
&(potewald[0][0][0]));
H5Dclose(h_data);
H5Sclose(h_space);
H5Fclose(h_file);
#endif
......@@ -207,6 +229,7 @@ void gravity_exact_force_ewald_init(double boxSize) {
fewald_x[i][j][k] *= boxSize_inv2;
fewald_y[i][j][k] *= boxSize_inv2;
fewald_z[i][j][k] *= boxSize_inv2;
potewald[i][j][k] *= boxSize_inv2;
}
}
}
......@@ -229,11 +252,12 @@ void gravity_exact_force_ewald_init(double boxSize) {
* @param rx x-coordinate of distance vector.
* @param ry y-coordinate of distance vector.
* @param rz z-coordinate of distance vector.
* @param corr (return) The Ewald correction.
* @param corr_f (return) The Ewald correction for the force.
* @param corr_p (return) The Ewald correction for the potential.
*/
__attribute__((always_inline)) INLINE static void
gravity_exact_force_ewald_evaluate(double rx, double ry, double rz,
double corr[3]) {
double corr_f[3], double *corr_p) {
const double s_x = (rx < 0.) ? 1. : -1.;
const double s_y = (ry < 0.) ? 1. : -1.;
......@@ -258,40 +282,51 @@ gravity_exact_force_ewald_evaluate(double rx, double ry, double rz,
const double tz = 1. - dz;
/* Interpolation in X */
corr[0] = 0.;
corr[0] += fewald_x[i + 0][j + 0][k + 0] * tx * ty * tz;
corr[0] += fewald_x[i + 0][j + 0][k + 1] * tx * ty * dz;
corr[0] += fewald_x[i + 0][j + 1][k + 0] * tx * dy * tz;
corr[0] += fewald_x[i + 0][j + 1][k + 1] * tx * dy * dz;
corr[0] += fewald_x[i + 1][j + 0][k + 0] * dx * ty * tz;
corr[0] += fewald_x[i + 1][j + 0][k + 1] * dx * ty * dz;
corr[0] += fewald_x[i + 1][j + 1][k + 0] * dx * dy * tz;
corr[0] += fewald_x[i + 1][j + 1][k + 1] * dx * dy * dz;
corr[0] *= s_x;
corr_f[0] = 0.;
corr_f[0] += fewald_x[i + 0][j + 0][k + 0] * tx * ty * tz;
corr_f[0] += fewald_x[i + 0][j + 0][k + 1] * tx * ty * dz;
corr_f[0] += fewald_x[i + 0][j + 1][k + 0] * tx * dy * tz;
corr_f[0] += fewald_x[i + 0][j + 1][k + 1] * tx * dy * dz;
corr_f[0] += fewald_x[i + 1][j + 0][k + 0] * dx * ty * tz;
corr_f[0] += fewald_x[i + 1][j + 0][k + 1] * dx * ty * dz;
corr_f[0] += fewald_x[i + 1][j + 1][k + 0] * dx * dy * tz;
corr_f[0] += fewald_x[i + 1][j + 1][k + 1] * dx * dy * dz;
corr_f[0] *= s_x;
/* Interpolation in Y */
corr[1] = 0.;
corr[1] += fewald_y[i + 0][j + 0][k + 0] * tx * ty * tz;
corr[1] += fewald_y[i + 0][j + 0][k + 1] * tx * ty * dz;
corr[1] += fewald_y[i + 0][j + 1][k + 0] * tx * dy * tz;
corr[1] += fewald_y[i + 0][j + 1][k + 1] * tx * dy * dz;
corr[1] += fewald_y[i + 1][j + 0][k + 0] * dx * ty * tz;
corr[1] += fewald_y[i + 1][j + 0][k + 1] * dx * ty * dz;
corr[1] += fewald_y[i + 1][j + 1][k + 0] * dx * dy * tz;
corr[1] += fewald_y[i + 1][j + 1][k + 1] * dx * dy * dz;
corr[1] *= s_y;
corr_f[1] = 0.;
corr_f[1] += fewald_y[i + 0][j + 0][k + 0] * tx * ty * tz;
corr_f[1] += fewald_y[i + 0][j + 0][k + 1] * tx * ty * dz;
corr_f[1] += fewald_y[i + 0][j + 1][k + 0] * tx * dy * tz;
corr_f[1] += fewald_y[i + 0][j + 1][k + 1] * tx * dy * dz;
corr_f[1] += fewald_y[i + 1][j + 0][k + 0] * dx * ty * tz;
corr_f[1] += fewald_y[i + 1][j + 0][k + 1] * dx * ty * dz;
corr_f[1] += fewald_y[i + 1][j + 1][k + 0] * dx * dy * tz;
corr_f[1] += fewald_y[i + 1][j + 1][k + 1] * dx * dy * dz;
corr_f[1] *= s_y;
/* Interpolation in Z */
corr[2] = 0.;
corr[2] += fewald_z[i + 0][j + 0][k + 0] * tx * ty * tz;
corr[2] += fewald_z[i + 0][j + 0][k + 1] * tx * ty * dz;
corr[2] += fewald_z[i + 0][j + 1][k + 0] * tx * dy * tz;
corr[2] += fewald_z[i + 0][j + 1][k + 1] * tx * dy * dz;
corr[2] += fewald_z[i + 1][j + 0][k + 0] * dx * ty * tz;
corr[2] += fewald_z[i + 1][j + 0][k + 1] * dx * ty * dz;
corr[2] += fewald_z[i + 1][j + 1][k + 0] * dx * dy * tz;
corr[2] += fewald_z[i + 1][j + 1][k + 1] * dx * dy * dz;
corr[2] *= s_z;
corr_f[2] = 0.;
corr_f[2] += fewald_z[i + 0][j + 0][k + 0] * tx * ty * tz;
corr_f[2] += fewald_z[i + 0][j + 0][k + 1] * tx * ty * dz;
corr_f[2] += fewald_z[i + 0][j + 1][k + 0] * tx * dy * tz;
corr_f[2] += fewald_z[i + 0][j + 1][k + 1] * tx * dy * dz;
corr_f[2] += fewald_z[i + 1][j + 0][k + 0] * dx * ty * tz;
corr_f[2] += fewald_z[i + 1][j + 0][k + 1] * dx * ty * dz;
corr_f[2] += fewald_z[i + 1][j + 1][k + 0] * dx * dy * tz;
corr_f[2] += fewald_z[i + 1][j + 1][k + 1] * dx * dy * dz;
corr_f[2] *= s_z;
/* Interpolation for potential */
*corr_p = 0.;
*corr_p += potewald[i + 0][j + 0][k + 0] * tx * ty * tz;
*corr_p += potewald[i + 0][j + 0][k + 1] * tx * ty * dz;
*corr_p += potewald[i + 0][j + 1][k + 0] * tx * dy * tz;
*corr_p += potewald[i + 0][j + 1][k + 1] * tx * dy * dz;
*corr_p += potewald[i + 1][j + 0][k + 0] * dx * ty * tz;
*corr_p += potewald[i + 1][j + 0][k + 1] * dx * ty * dz;
*corr_p += potewald[i + 1][j + 1][k + 0] * dx * dy * tz;
*corr_p += potewald[i + 1][j + 1][k + 1] * dx * dy * dz;
}
#endif
......@@ -383,6 +418,7 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
/* Be ready for the calculation */
double a_grav[3] = {0., 0., 0.};
double pot = 0.;
/* Interact it with all other particles in the space.*/
for (int j = 0; j < (int)s->nr_gparts; ++j) {
......@@ -408,37 +444,42 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
const double r_inv = 1. / sqrt(r2);
const double r = r2 * r_inv;
const double mj = gpj->mass;
double f;
double f, phi;
if (r >= hi) {
/* Get Newtonian gravity */
f = mj * r_inv * r_inv * r_inv;
phi = -mj * r_inv;
} else {
const double ui = r * hi_inv;
double W;
double Wf, Wp;
kernel_grav_eval_double(ui, &W);
kernel_grav_eval_force_double(ui, &Wf);
kernel_grav_eval_pot_double(ui, &Wp);
/* Get softened gravity */
f = mj * hi_inv3 * W;
f = mj * hi_inv3 * Wf;
phi = mj * hi_inv * Wp;
}
a_grav[0] += f * dx;
a_grav[1] += f * dy;
a_grav[2] += f * dz;
pot += phi;
/* Apply Ewald correction for periodic BC */
if (periodic && r > 1e-5 * hi) {
double corr[3];
gravity_exact_force_ewald_evaluate(dx, dy, dz, corr);
double corr_f[3], corr_pot;
gravity_exact_force_ewald_evaluate(dx, dy, dz, corr_f, &corr_pot);
a_grav[0] += mj * corr[0];
a_grav[1] += mj * corr[1];
a_grav[2] += mj * corr[2];
a_grav[0] += mj * corr_f[0];
a_grav[1] += mj * corr_f[1];
a_grav[2] += mj * corr_f[2];
pot += mj * corr_pot;
}
}
......@@ -446,6 +487,7 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
gpi->a_grav_exact[0] = a_grav[0] * const_G;
gpi->a_grav_exact[1] = a_grav[1] * const_G;
gpi->a_grav_exact[2] = a_grav[2] * const_G;
gpi->potential_exact = pot * const_G;
counter++;
}
......@@ -529,8 +571,9 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
fprintf(file_swift, "# theta= %16.8e\n", e->gravity_properties->theta_crit);
fprintf(file_swift, "# Git Branch: %s\n", git_branch());
fprintf(file_swift, "# Git Revision: %s\n", git_revision());
fprintf(file_swift, "# %16s %16s %16s %16s %16s %16s %16s\n", "id", "pos[0]",
"pos[1]", "pos[2]", "a_swift[0]", "a_swift[1]", "a_swift[2]");
fprintf(file_swift, "# %16s %16s %16s %16s %16s %16s %16s %16s\n", "id",
"pos[0]", "pos[1]", "pos[2]", "a_swift[0]", "a_swift[1]",
"a_swift[2]", "potential");
/* Output particle SWIFT accelerations */
for (size_t i = 0; i < s->nr_gparts; ++i) {
......@@ -541,9 +584,10 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
gpart_is_starting(gpi, e)) {
fprintf(file_swift, "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e \n",
fprintf(file_swift,
"%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",
gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2]);
gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2], gpi->potential);
}
}
......@@ -569,9 +613,9 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
fprintf(file_exact, "# periodic= %d\n", s->periodic);
fprintf(file_exact, "# Git Branch: %s\n", git_branch());
fprintf(file_exact, "# Git Revision: %s\n", git_revision());
fprintf(file_exact, "# %16s %16s %16s %16s %16s %16s %16s\n", "id",
fprintf(file_exact, "# %16s %16s %16s %16s %16s %16s %16s %16s\n", "id",
"pos[0]", "pos[1]", "pos[2]", "a_exact[0]", "a_exact[1]",
"a_exact[2]");
"a_exact[2]", "potential");
/* Output particle exact accelerations */
for (size_t i = 0; i < s->nr_gparts; ++i) {
......@@ -582,10 +626,11 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
gpart_is_starting(gpi, e)) {
fprintf(
file_exact, "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e \n",
gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2]);
fprintf(file_exact,
"%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",
gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
gpi->a_grav_exact[0], gpi->a_grav_exact[1],
gpi->a_grav_exact[2], gpi->potential_exact);
}
}
......
src/gravity/Default/gravity.h
View file @ df794475
......@@ -48,16 +48,29 @@ __attribute__((always_inline)) INLINE static float gravity_get_softening(
}
/**
* @brief Returns the potential of a particle
* @brief Returns the comoving potential of a particle
*
* @param gp The particle of interest
*/
__attribute__((always_inline)) INLINE static float gravity_get_potential(
const struct gpart* restrict gp) {
__attribute__((always_inline)) INLINE static float
gravity_get_comoving_potential(const struct gpart* restrict gp) {
return gp->potential;
}
/**
* @brief Returns the physical potential of a particle
*
* @param gp The particle of interest.
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static float
gravity_get_physical_potential(const struct gpart* restrict gp,
const struct cosmology* cosmo) {
return gp->potential * cosmo->a_inv;
}
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
......@@ -136,6 +149,7 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
gp->a_grav[0] *= const_G;
gp->a_grav[1] *= const_G;
gp->a_grav[2] *= const_G;
gp->potential *= const_G;
}
/**
......
src/gravity/Default/gravity_iact.h
View file @ df794475
......@@ -38,9 +38,11 @@
* @param h_inv3 Cube of the inverse of the softening length.
* @param mass Mass of the point-mass.
* @param f_ij (return) The force intensity.
* @param pot_ij (return) The potential.
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij) {
float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij,
float *pot_ij) {
/* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2);
......@@ -50,17 +52,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
/* Get Newtonian gravity */
*f_ij = mass * r_inv * r_inv * r_inv;
*pot_ij = -mass * r_inv;
} else {
const float r = r2 * r_inv;
const float ui = r * h_inv;
float W_ij;
kernel_grav_eval(ui, &W_ij);
float W_f_ij, W_pot_ij;
kernel_grav_force_eval(ui, &W_f_ij);
kernel_grav_pot_eval(ui, &W_pot_ij);
/* Get softened gravity */
*f_ij = mass * h_inv3 * W_ij;
*f_ij = mass * h_inv3 * W_f_ij;
*pot_ij = mass * h_inv * W_pot_ij;
}
}
......@@ -78,10 +83,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
* @param mass Mass of the point-mass.
* @param rlr_inv Inverse of the mesh smoothing scale.
* @param f_ij (return) The force intensity.
* @param pot_ij (return) The potential.
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
float r2, float h2, float h_inv, float h_inv3, float mass, float rlr_inv,
float *f_ij) {
float *f_ij, float *pot_ij) {
/* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2);
......@@ -92,23 +98,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
/* Get Newtonian gravity */
*f_ij = mass * r_inv * r_inv * r_inv;
*pot_ij = -mass * r_inv;
} else {
const float ui = r * h_inv;
float W_ij;
float W_f_ij, W_pot_ij;
kernel_grav_eval(ui, &W_ij);
kernel_grav_force_eval(ui, &W_f_ij);
kernel_grav_pot_eval(ui, &W_pot_ij);
/* Get softened gravity */
*f_ij = mass * h_inv3 * W_ij;</