Commit 61dcf239 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Make the gravity kernel testing script closer to the Gadget code.

parent 23741d2d
......@@ -31,26 +31,31 @@
/**
* @brief The Gadget-2 gravity kernel function
*
* Taken from Gadget-2.0.7's forcetree.c lines 2755-2800
*
* @param r The distance between particles
* @param h The cut-off distance of the kernel
*/
float gadget(float r, float h) {
float fac;
const float r2 = r * r;
if (r >= h)
fac = 1.0f / (r2 * r);
else {
const float h_inv = 1. / h;
const float h_inv3 = h_inv * h_inv * h_inv;
const float u = r * h_inv;
float gadget(float r, float epsilon) {
const float h = epsilon;
const float h_inv = 1.f / h;
const float u = r * h_inv;
if (u >= 1) {
const float r_inv = 1. / r;
return r_inv * r_inv * r_inv;
} else {
if (u < 0.5)
fac = h_inv3 * (10.666666666667 + u * u * (32.0 * u - 38.4));
return h_inv * h_inv * h_inv *
(10.666666666667 + u * u * (32.0 * u - 38.4));
else
fac =
h_inv3 * (21.333333333333 - 48.0 * u + 38.4 * u * u -
10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
return h_inv * h_inv * h_inv *
(21.333333333333 - 48.0 * u + 38.4 * u * u -
10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
}
return fac;
}
int main() {
......@@ -61,20 +66,22 @@ int main() {
for (int k = 1; k < numPoints; ++k) {
const float r = (r_max * k) / numPoints;
const float u = r / h;
const float gadget_w = gadget(r, h);
const float h_inv = 1.f / h;
const float h_inv3 = h_inv * h_inv * h_inv;
const float u = r * h_inv;
float swift_w;
if (u < 1.) {
kernel_grav_eval(u, &swift_w);
swift_w *= (1 / (h * h * h));
} else {
if (r >= h) {
swift_w = 1 / (r * r * r);
} else {
kernel_grav_eval(u, &swift_w);
swift_w *= h_inv3;
}
if (fabsf(gadget_w - swift_w) > 2e-7) {
if (fabsf(gadget_w - swift_w) > 1e-5 * fabsf(gadget_w)) {
printf("%2d: r= %f h= %f u= %f Wg(r,h)= %f Ws(r,h)= %f\n", k, r, h, u,
gadget_w, swift_w);
......@@ -87,28 +94,30 @@ int main() {
printf("\nAll values are consistent\n");
/* Now test the long range function */
const float a_smooth = 4.5f;
/* const float a_smooth = 4.5f; */
for (int k = 1; k < numPoints; ++k) {
/* for (int k = 1; k < numPoints; ++k) { */
const float r = (r_max * k) / numPoints;
/* const float r = (r_max * k) / numPoints; */
const float u = r / a_smooth;
/* const float u = r / a_smooth; */
float swift_w;
kernel_long_grav_eval(u, &swift_w);
/* float swift_w; */
/* kernel_long_grav_eval(u, &swift_w); */
float gadget_w = erfc(u / 2) + u * exp(-u * u / 4) / sqrt(M_PI);
/* float gadget_w = erfcf(u / 2) + u * expf(-u * u / 4) / sqrtf(M_PI); */
if (fabsf(gadget_w - swift_w) > 2e-7) {
/* if (fabsf(gadget_w - swift_w) > 1e-4 * fabsf(gadget_w)) { */
printf("%2d: r= %f r_lr= %f u= %f Ws(r)= %f Wg(r)= %f\n", k, r, a_smooth,
u, swift_w, gadget_w);
/* printf("%2d: r= %f r_lr= %f u= %f Ws(r)= %f Wg(r)= %f\n", k, r,
* a_smooth, */
/* u, swift_w, gadget_w); */
printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w);
return 1;
}
}
/* printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w);
*/
/* return 1; */
/* } */
/* } */
return 0;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment